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Introduction 


The  hypothalamic-pituitary-adrenal  (HPA)  axis  controls  the  body’s  “fight  or  flight”  response  through  a  series 
of  endocrine  and  immune  signals  directed  at  ensuring  immediate  survival  and  later  re-establishing 
homeostasis.  Changes  in  the  tone  of  this  response  have  been  observed  in  veterans  with  Gulf  War  Illness 
(GWI).  Studies  report  abnormal  cell  proliferation,  impaired  function  and  persistent  oxidative  stress  in 
circulating  immune  cells  of  patients.  Similarly  dysregulation  of  the  HPA  axis  includes  hypersensitivity  in 
cytokine  feedback  as  well  as  suppression  of  cortisol  and  neurotransmitters  responsible  for  mediating  innate 
and  adaptive  immunity.  This  is  further  complicated  by  the  impact  on  the  HPA  axis  of  a  myriad  of  regulatory 
interactions  both  within  and  between  i)  the  immune  system  and  ii)  the  sex-hormone  axis,  the  hypothalamic- 
pituitary-gonadal  (HPG)  axis. 

We  proposed  that  severe  physical  or  psychological  insult  to  the  endocrine  and  immune  systems  can 
displace  these  from  a  normal  regulatory  equilibrium  into  a  compromised  stable  state.  This  state  is 
characterized  by  a  self-perpetuating  inflammatory  response  that  involves  regulatory  imbalance  between  the 
HPA,  HPG  and  immune  axes.  To  explore  the  validity  of  this  hypothesis  our  objective  was  to  create 
comprehensive  engineering  models  of  endocrine-immune  interaction  dynamics  in  order  to  identify  (i) 
theoretical  failure  modes  of  the  endocrine-immune  interplay  that  align  with  GWI,  and  (ii)  promising  treatment 
strategies  that  exploit  the  naturally  occurring  stable  points  of  these  systems. 

Body. 

At  the  time  of  our  last  update  we  had  completed  a  re-assessment  of  the  modeling  approach  and 
successfully  identified,  refined  and  deployed  a  discrete  modeling  paradigm  enabling  us  to  circumvent  the 
significant  gaps  in  the  required  parameter  estimates  exist  in  the  literature.  This  new  approach  was 
described  in  greater  detail  in  the  previous  report  (September,  2012)  (Task  1)  and  consists  in  an  extension 
of  the  discrete  logical  network  methodology  proposed  originally  by  Thomas  et  al.  [1 ,2]  and  developed  further 
by  Mendoza  and  Xenarios  [3],  Importantly,  this  approach  supports  the  seamless  integration  of  kinetic 
information  wherever  available,  be  it  simple  sequential  precedence,  relative  time  scale  or  detailed 
dynamics. 

We  have  now  completed  the  original  scope  of  Task  1  through  Task  5.  However  as  a  result  of  improvements 
in  computational  efficiency  we  have  been  able  to  extend  the  basic  models  beyond  what  was  originally 
anticipated.  With  the  concerted  efforts  of  Dr.  Craddock,  new  programming  staff  Mark  Rice  and  Ryan  del 
Rosario  and  research  interns  Simar  Singh  and  Lundy  McKibbin  we  have  continued  to:  (1)  design  and 
deploy  a  distributed  version  of  the  computing  code,  (2)  increase  the  granularity  of  the  multi-axis  model  and 
the  implement  a  statistical  scheme  for  model  validation,  (3)  significantly  extend  the  scope  and  increase 
fidelity  of  the  detailed  immune  model,  and  (4)  develop  a  prototype  model  of  neuro-inflammation,  and  (5) 
designed  and  deployed  a  first  prototype  of  the  treatment  optimization  scheme. 

1.  Continued  algorithm  development  and  speed-up.  (Task  1).  The  core  concept  of  the  approach  we 
have  used  is  connectivity.  Key  biological  regulatory  processes  have  been  translated  into  a  set  of  discrete 
logic  circuits.  Analysis  of  these  networks  makes  it  possible  to  identify  the  number  and  type  (e.g.  oscillatory, 
etc...)  of  resting  states  as  well  as  their  molecular  and  cellular  profile  without  detailed  knowledge  of  response 
dynamics.  Early  implementations  of  this  analysis  were  made  in  a  high-level  rapid-prototyping  environment 
(Python)  facilitating  development  but  severely  limiting  computational  performance.  As  mentioned  previously, 
under  this  discrete  formalism  the  number  of  model  variables  determines  the  total  number  of  system-wide 
states  such  that  a  model  of  N  state  variables  possesses  3N  states.  As  a  result  the  number  of  total  system- 
wide  states  increases  rapidly  as  new  state  variable  elements  are  added.  Initially  these  calculations  were 
encoded  into  a  rapid-prototyping  Python  script  that  was  used  to  search  the  above-mentioned  network  for 
stable  equilibrium  states  (Version  0).  Within  a  24-hour  threshold  time  (86,400  seconds)  this  version  is 
capable  of  analyzing  up  to  14  variables  (4,782,969  states)  with  a  memory  usage  in  the  range  of  gigabytes 
(GB)  (Figure  1).  High-performance  computing  staff,  programmers  Rice  and  del  Rosario,  re-engineered  the 
search  algorithm  and  its  implementation  in  several  stages.  First,  the  algorithm  was  directly  re-coded  in  the 
C  programming  language  (Version  1 )  as  it  is  both  memory-efficient  and  approximately  30  times  faster  than 
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Python.  This  increased  performance  enabled  analysis  of  a  17-variable  model  (129,140,163  possible 
discrete  states)  in  24  hours  at  a  memory  usage  in  the  megabyte  (MB)  range.  Next,  the  serial  algorithm  was 
re-engineered  through  the  introduction  of  ternary  data  structures  to  efficiently  optimize  memory  usage  and 
run  time,  and  prepare  the  algorithm  for  parallel  implementation  (Version  2).  Here,  in  the  24-hour  threshold 
time,  a  19-variable  model  (1,162,261,467  states)  was  analyzed  at  MB  memory  usage.  Thirdly,  the  algorithm 
was  implemented  with  parallel  tasking  (Version  3),  using  multiple  levels  of  parallel  threading  (mO  to  m4). 
Here  we  have  successfully  run  a  model  with  23  variables  (94,143,178,827  states)  within  a  day  using  only 
MB’s  of  memory.  Finally,  we  have  parallelized  the  code  further  with  a  supervisory  layer  based  on  message 
passing  interface  (MPI)  (Version  4)  to  make  full  use  of  the  high-performance  computing  resources  on  the 
University  of  Miami’s  Pegasus  cluster.  While  performance  measures  are  still  being  evaluated  we  have 
successfully  analyzed  a  25  variable  circuit  model  (847,288,609,443  discrete  states).  Continued 
improvements  to  computational  efficiency  are  ongoing. 


2.  Continued  refinement  of  an  integrated  model  of  HPA-HPG-immune  interaction  ( Task  3,  Task  5). 

We  had  previously  extended  our  early  model  of  HPA  axis  dynamics  [4]  by  including  feed-forward  and 
feedback  interactions  with  sex  hormone  regulation  and  immune  response.  A  circuit  model  had  been 
constructed  that  linked  state  variables  across  the  HPA  axis  with  hypothalamic-pituitary-gonadal  (HPG) 
function  in  both  men  and  women,  as  well  as  a  coarse-grained  mode  of  the  immune  system  consisting  of 
innate  (MR)  and  adaptive  (AIR)  immune  components.  A  critical  review  of  this  model  prompted  us  to:  (i)  re¬ 
assess  the  coarse-graining  of  the  immune  components  (MR  and  AIR  aggregate  nodes),  increasing  the  level 
of  detail  to  improve  fidelity,  and  (ii)  define  and  implement  an  alternate  validation  measure. 

•  A  modified  multi-axis  model.  In  a  first  coarse  iteration  of  the  model,  immune  function  was  described 
simply  in  terms  cytokine  activity  of  the  innate  (HR)  and  adaptive  (AIR)  immune  responses.  Here  the 
aggregation  of  all  adaptive  immune  response  into  the  AIR  node  lacked  the  complexity  needed  to  capture 
shifts  between  Thl  and  Th2  activity.  Further  resolution  was  added  to  the  immune  model  by  separating 
the  AIR  into  Thl  and  Th2  activity  and  by  adding  both  cell  population  activity,  as  well  as  cytokine 
signaling  separately.  In  this  modified  immune  module  innate  immune  cells  (ICells)  produce  cytokines 
that  regulate  the  innate  immune  response  (HR)  including  interleukin  (IL)  -1,  IL-6,  IL-8,  IL-12,  IL-15,  IL-23 
and  tumor  necrosis  factor  alpha  (TNF-a).  These  HR  signals  serve  to  prime  helper  T  cells  towards  a  Thl 
type  adaptive  immune  response  (T 1  Cell),  producing  Thl  pro-inflammatory  cytokines  (TICyt)  including 
IL-2,  interferon-gamma  (IFN-y),  and  tumor  necrosis  factor  beta  (TNF-P).  This  further  activates  ICells, 
while  suppressing  the  Th2  adaptive  immune  response  (T2Cell).  The  T2Cell  node  promotes  the 
production  of  the  Th2  anti-inflammatory  cytokines  (T2Cyt)  IL-4,  IL-5,  IL-10  and  IL-13,  which  serve  to 
inhibit  the  activity  of  T1  Cell  and  ICells.  Interaction  with  the  HPA  axis  is  mediated  by  CORT  suppressing 
ICell  and  Th  1  Cell  activity,  while  HR  and  ThICyt  signals  stimulate  the  HPA.  Additionally,  new  interactions 
between  the  immune  and  HPG  axis  were  included  where  ThICyt  signals  suppress  GnRH  and  LH/FSH 
release  and  the  dimorphic  response  of  sex  hormones  TEST/  EST  serve  to  induce  Th1/Th2  activity. 

•  A  probabilistic  measure  of  alignment  with  experimental  data.  Alignment  of  model  predictions  with 
experimental  data  were  previously  assessed  on  the  basis  of  discrete  Hamming  distance  and  visualized 
with  a  Sammon  projection  of  the  latter.  In  order  to  provide  a  more  continuous  measure  of  similarity  or 
dissimilarity  we  have  adopted  a  probabilistic  measure  proposed  by  Brown  [5],  Here,  we  calculate  the 
significance  of  alignment  between  experimental  data  and  a  given  state  predicted  by  the  model  using  a 
meta-analysis  technique  that  combines  non-independent  test  statistics.  Null  probability  p-values  for 
individual  variables  are  calculated  using  two-sample  t-tests  between  ill  subjects  and  healthy  controls. 

To  give  the  probability  of  obtaining  the  model  value  by  chance  ‘right-handed’  one-tailed  tests  are  used 
when  the  model  predicts  a  high  state,  ‘left-handed’  tests  when  predictions  are  low,  and  two-tailed  tests 
when  the  prediction  is  a  nominal  value.  These  non-independent  statistics  are  then  combined  into  a  chi- 
squared  test  statistic,  which  is  scaled  to  T  =  T0/c  with  2N/c  degrees  of  freedom,  where  c  =  ct2/4N.  This 
statistic  is  then  used  in  the  scaled  chi-squared  distribution  to  determine  the  overall  probability  of 
obtaining  the  alignment  by  chance.  The  advantage  of  this  method  is  that  it  accommodates  for  the 
dependence  between  variables,  allows  for  a  statement  of  confidence  on  alignment  for  each  individual 
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model  predicted  state,  and  does  not  depend  on  the  number  of  measureable  markers  allowing  for  direct 
comparison  across  models. 

Details  of  this  analysis  and  the  final  multi-axis  model  are  described  in  Appendix  A  in  manuscript  recently 
submitted  to  PLoS  One  [6],  In  brief,  co-regulation  of  the  HPA,  HPG  and  immune  systems  has  been 
described  as  a  revised  circuit  model  consisting  of  14  state  variables  where  each  variable  can  assume  one 
of  three  discrete  states  at  any  point  in  time:  -1  (inhibited),  0  (nominal)  and  +1  (elevated).  In  this  model  the 
HPA  axis  continues  to  be  described  in  terms  of  corticotropin-releasing  hormone  (CRH),  adrenocorticotropic 
hormone  (ACTH),  cortisol  (CORT)  and  cytostolic  glucocorticoid  receptors  (GRs),  which  unlike  membrane 
bound  receptors,  dimerize  (GRD)  (Figure  2A-B).  HPG  function  is  again  described  by  the  levels  of 
gonadotropin-releasing  hormone  (GnRH),  of  luteinizing  homone  (LH)  as  well  as  testosterone  (TEST)  in 
males  (Figure  2  C)  and  estradiol  (EST)  in  females  (Figure  2  D-G).  As  before,  the  effects  of  gender  merit 
special  attention  as  testosterone  (TEST)  exhibits  an  inhibitory  effect  on  the  HPA  axis  while  estrogen  (EST) 
and  progesterone  can  stimulate  or  suppress  HPA  activity  depending  on  the  phase  of  menstrual  cycle. 
Theses  components  are  integrated  with  the  revised  immune  model  described  above. 

Results  of  simulations  conducted  on  the  refined  model  can  be  summarized  as  follows: 

•  Male  subjects.  Inclusion  of  basic  immune  function  and  sex  hormone  regulation  by  the  HPG  axis  with 
HPA  function  (HPA-GR-lmmune-HPG  model)  in  male  subjects  (Figure  2C)  resulted  in  the  emergence  of 
5  stable  equilibrium  states.  Once  again  the  first  state  was  that  of  normal  health  (SSO).  Low  levels  of 
ACTH  and  elevated  expression  of  the  glucocorticoid  receptors  GR  and  GRD  characterized  the  second 
equilibrium  state  (SSI).  The  third  stable  state  (SS2)  exhibited  supprssed  innate  and  Thl  immune 
responses  (low  ICell,  HR,  T 1  Cell,  and  TICyt),  with  increased  Th2  activity  (high  T2Cell  and  T2Cyt).  The 
fourth  state  (SS3)  presented  low  ACTH,  suppressed  innate  and  Thl  immune  activity  (low  ICell,  HR, 

T 1  Cell  and  TICyt),  and  elevated  Th2  and  glucocorticoid  receptor  activity  (high  GRD,  GR,  T2Cell  and 
T2Cyt).  The  final  state  (SS4)  displayed  hypercortisolism,  suppressed  HPG  activity  and  a  shift  towards 
the  Thl  immune  response  (low  T2Cell,  T2Cyt,  GnRH,  LH/FSH  and  TEST/EST,  and  high  CORT,  GRD, 
GR,  TICyt  and  TICell). 

•  Female  subjects.  In  the  specific  case  of  positive  feedback  along  the  HPG  axis  and  suppressive 
interaction  with  the  HPA  axis,  the  HPA-GR-lmmune-HPG  model  for  female  subjects  (Figure  2F) 
supported  1 1  steady  states.  In  addition  to  5  states  equivalent  to  those  obtained  for  the  male  subjects, 
we  found  new  steady  states  that  corresponded  to  suppressed  HPA  axis  and  innate  immune  response 
(low  CRH,  ACTH,  CORT,  ICell  and  HR),  while  the  HPG  and  anti-inflammatory  response  were  elevated 
(high  T2Cell,  T2Cyt,  GnRH,  LH/FSH  and  EST).  This  combination  occurred  at  each  of  the  three  low, 
nominal  and  high  values  for  glucocorticoid  receptor  activity  (GR/GRD)  (SS5,  SS6  and  SS7, 
respectively).  The  final  three  additional  states  all  supported  suppressed  HPA  (CRH,  ACTH,  and  CORT) 
and  TICell  activity,  with  elevated  HPG  activity  (GnRH,  LH/FSH  and  EST).  These  were  again 
differentiated  by  their  glucocorticoid  receptor  levels  (GR/GRD  at  low  (SS8),  nominal  (SS9)  and  high 
(SS10)  values).  Note  that  a  stable  steady  state  characterized  by  low  cortisol  levels  was  found  only  for 
female  subjects. 

•  Alignment  with  experimental  data.  To  validate  these  results  the  predicted  steady  states  were  first 
compared  to  steroid  and  cytokine  levels  recorded  in  male  Gulf  War  veterans  with  GWI  and  healthy 
veterans  (HCs)  as  part  of  a  sister  study  [7].  As  experimental  measures  for  ACTH,  GR,  GRD,  and 
immune  cells  populations  were  not  available,  certain  steady  states  could  not  be  distinguished  and 
validated  separately  from  one  another.  Comparison  to  the  nominal  states  (SS0/SS1)  showed  poor 
alignment,  with  a  null  probability  of  p=0.82,  suggesting  that  the  GWI  profile  cannot  be  considered  the 
same  as  normal  behavior.  The  predicted  states  presenting  a  shift  towards  Th2  immune  activation 
(SS2/SS3)  showed  improved  alignment  with  a  significance  of  p=0.38.  However  the  final  state  (SS4) 
displaying  hyper-cortisolism,  low  TEST  and  a  shift  towards  Thl  immune  activation  yielded  the  best 
alignment  with  a  null  probability  of  p=0.30,  supporting  the  notion  of  a  more  classical  Thl  auto-immune 
signature  with  a  concurrent  (and  perhaps  stabilizing)  endocrine  component  in  GWI. 

As  a  much  greater  proportion  of  women  than  men  are  affected  by  CFS,  we  compared  the  predicted 
steady  states  identified  with  the  female  model  to  experimental  data  collected  under  two  compatible 
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studies  [8,9],  Alignment  with  the  baseline  nominal  setting  in  measureable  variables  (SSO/SS1)  was  poor, 
p=0.83,  reinforcing  that  CFS  is  distinctly  different  from  normal  regulatory  behavior.  The  Th2-shifted 
immune  profiles  predicted  by  the  model  (SS2/SS3)  showed  a  significant  alignment  with  the  measured 
signature  (p=0.04),  suggesting  that  Th2  activation  in  CFS  may  at  least  in  part  be  supported  by 
homeostatic  drive.  This  emphasized  further  by  a  low  degree  of  alignment  with  the  Thl  immune  activated 
state,  with  hypercortisolism,  and  low  EST  (p=0.28).  Improved  alignment  was  found  with  states  with  a 
shift  towards  Th2,  coupled  with  hypocortisolism,  and  high  EST  (SS5/SS6/SS7)  (p=0.02).  States 
presenting  with  only  hypocortisolism  and  high  EST,  and  no  immune  activation  (SS8/SS9/SS10)  aligned 
very  weakly  with  the  measured  profile  (p=0.60),  suggesting  that  hypocortisolism,  increased  EST,  and 
Th2  activation  in  combination  are  key  CFS  profile  features  that  might  owe  at  least  part  of  their 
persistence  to  basic  homeostatic  control. 


3.  Refinement  of  detailed  immune  circuitry  including  Thl  7  and  neurotransmission  (Task  2,  3  and  5). 
Based  on  the  work  of  Folcik  et  al.  (2007,  201 1 )  [1 0, 1 1  ]  and  an  extensive  review  of  recent  literature,  we 
constructed  an  initial  wiring  diagram  describing  cytokine  signaling  between  immune  cell  populations  [12] 

(see  excerpts  of  manuscript  in  preparation;  Appendix  B  in  2012  annual  report).  We  have  now  extended  this 
first  detailed  model  of  immune  signaling  at  two  levels  of  granularity: 

•  Addition  of  Tree  and  Thl  7  components  to  aggregate  model  of  cytokine  signaling.  The  specific 

cytokines  supported  in  the  initial  model  included  interleukin  (IL)-1,  IL-2,  IL-4,  IL-5,  IL-6,  IL-8,  IL-10,  IL-12, 
IL-13,  IL-23,  IL-27,  interferon  (IFN)-y,  and  tumor  necrosis  factor  (TNF)-a.  In  order  to  improve 
computational  efficiency,  cytokines  were  grouped  according  to  their  dominant  action  into  either  a 
monokine  (MK)  or  cytokine  (CK)  group  (Folcik  et  al.,  201 1  )[1 1],  We  have  now  extended  the  model  to 
include  the  actions  of  IL-17,  21  and  23  as  well  as  TGF-p.  The  extended  model  also  includes  T  regulatory 
(Treg),  Th17  and  activated  Th17(23)  cell  populations  as  well  as  IgA  and  IgG  antibody  classes.  This 
updated  version  of  the  cytokine  signaling  network  includes  as  before  the  effects  of  stress  and  sex 
hormones  for  male  subjects  only  at  this  time  (Figure  3).  Results  of  our  stability  analysis  on  this  second- 
generation  model  that  can  be  summarized  as  follows: 

o  Stable  immune  response  modes  in  male  subjects.  Application  of  this  discrete  dynamical 
analysis  to  the  detailed  endocrine-immune  network  yielded  three  predicted  stable  steady  states. 
As  always,  the  first  state  was  that  of  normal  health  (SS0).  The  second  stable  state  (SSI ) 
presented  with  low  anti-inflammatory  cytokines  (CK2),  low  testosterone,  and  suppressed  NK  cell 
activity,  Thl  and  Th2  immune  cell  activity,  accompanied  by  elevated  Thl  inflammatory  cytokines 
(CK1),  high  cortisol  levels,  and  increased  cytotoxic  T  lymphocyte  (CTL)  and  Treg  cell  activity. 
The  third  stable  attractor  (SS2)  also  displayed  low  NK  ceil  activity  and  testosterone  levels,  with 
elevated  cortisol,  CK1  and  CTL  activity.  However  this  state  presented  a  different  immune  profile 
characterized  by  low  TGF-p  levels,  elevated  monocyte  cytokines  (MK1,  MK2,  MK6,  MK21),  and 
increased  Th17  activity  (CK17,  Th17(23),  TH17b).  These  results  are  consistent  with  the  findings 
of  the  integrated  HPA-HPG-immune  model  discussed  above,  but  with  added  resolution  in  terms 
of  immune  function,  which  indicates  alternate  equilibria  defined  by  different  stable  levels  of 
regulatory  T  cell  activity,  or  Thl  7  immune  response. 

o  Alignment  with  experimental  data.  Once  again,  these  predicted  steady  states  were 

compared  with  experimental  data  used  in  Section  1  [7-9].  We  found  alignment  of  GWI  with  the 
healthy  reference  state  (SS0)  corresponded  to  a  null  probability  of  p=0.87,  indicating  very  poor 
alignment.  Improved  alignment  was  observed  with  SSI  (p=  0.30).  This  is  comparable  to  the 
results  found  for  the  high  CORT,  low  TEST,  elevated  Thl  response  state  using  the  integrated 
HPA-HPG-lmmune  model  described  previously  (section  2).  Further  improvement  was  found 
when  comparing  GWI  to  the  final  predicted  stable  state  (SS2)  (p=0.12),  suggesting  that  chronic 
Th17  activation,  hypercortisolism,  and  low  TEST  observed  in  this  illness  may  persist  in  part  as  a 
result  of  homeostatic  drive.  For  male  CFS  subjects  we  found  comparably  poor  alignment  with 
the  reference  baseline  state  SS0  (p=0.76).  However,  alignment  with  state  SSI  (p=0.16)  was 
dramatically  different  from  GWI,  emphasizing  the  distinct  nature  of  CFS.  Distinct  as  they  may 
be,  these  illnesses  nonetheless  share  some  common  components  that  our  group  has  begun  to 
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delineate  at  the  level  of  specific  pathways  [13],  Consistent  with  this  we  found  comparable 
alignment  of  CFS  and  GWI  with  SS2  (p=0.12)  albeit  for  slightly  different  reasons.  Thus,  while 
GWI  aligns  best  with  Th17  dysregulation,  CFS  alignment  suggests  slightly  different  imbalance  of 
Th17  and/or  Treg  response. 

•  High  fidelity  model  of  individual  cytokine  actions.  Improved  alignment  with  the  clinical  data  can  be 
accomplished  by  including  additional  key  interactions  in  the  model  regulatory  network,  and  by  increasing 
resolution  of  the  model  in  terms  of  the  state  variables  represented.  Key  interactions  located  outside  the 
immune  network  include  critical  neurotransmitters  linking  the  brain  and  central  nervous  system  with  the 
HPA  axis  and  the  immune  system.  The  neurotransmitters  norepinephrine  (NorEpi),  and  acetylcholine 
(ACh)  are  significant  regulators  of  cytokine  production,  and  therefore  immune  function.  Neuropeptide  Y 
(NPY)  is  also  key  component  of  the  stress  response,  and  its  subsequent  effects  on  the  immune  system. 
The  latter  has  now  been  shown  to  play  a  significant  role  in  CFS  [14].  These  messengers  have  now 
been  included  in  a  more  refined  model  of  the  immune  system  and  it’s  interface  with  neurotranmission 
(Figure  4).  Their  effect  on  the  immune  system  however  is  not  simple.  In  the  previous  immune  model 
several  cytokines  were  aggregated  into  groups.  NorEpi,  ACh  and  NPY  were  found  to  have  differing 
effects  on  the  production  of  cytokines  within  individual  groups.  To  accommodate  these  varied 
responses  several  of  the  aggregate  nodes  were  separated  into  their  individual  constituent  entities 
increasing  the  resolution,  and  complexity,  of  the  extended  immune  model.  The  overall  resulting  high 
fidelity  model  of  the  extended  immune  system,  including  HPA,  HPG  and  CNS  inputs  is  shown  in  Figure 
4.  A  first  analysis  has  shown  the  following: 

o  Preliminary  stable  immune  response  modes  in  male  subjects.  Discrete  logical  analysis  of  the 
preliminary  high  fidelity  extended  immune  model  produced  two  steady  states.  Normal  health 
characterized  the  first  state  (SSO),  while  the  second  state  (SSI)  presented  with  low  IL-1,  MK6, 
TEST  and  NK  cell  activity,  and  high  IL-12,  MK2,  CK1,  CK2,  CTL,  CORT  and  NPY.  Note  that 
several  of  the  cytokines  that  were  once  aggregated  (IL-1,  IL-8,  IL-12)  now  present  with  differing 
profiles.  These  yield  a  complicated  mixed  Thl  :Th2  profile  consistent  with  our  previous  analysis 
of  GWI.  Further  refinement  of  other  aggregate  nodes,  and  inclusion  of  the  Th17  axis  is  currently 
underway. 

o  Preliminary  alignment  with  experimental  data.  We  found  in  GWI  that  aligns  with  the  nominal 
steady  state  (SSO)  at  a  significance  level  p=0.85  again  indicating  poor  alignment  with  normal 
health.  Alignment  with  the  alternate  steady  state  (SSI )  however  was  much  more  significant  (p= 
0.07).  This  both  supports  a  notion  of  a  complex  stable  combined  Th1:Th2  response  in  this 
illness,  and  suggests  a  brain  component  in  its  perpetuation.  Further  analysis  with  the  refined 
model  is  being  conducted. 

Collectively  these  simulations  of  known  endocrine-immune  circuitry  support  the  existence  alternate 
homeostatic  regimes,  some  of  which  overlap  substantially  with  observed  immune  and  endocrine  status  in 
male  GWI  and  female  CFS  subjects.  Such  overlap  with  naturally  occurring  stable  regulatory  regimes  would 
certainly  be  consistent  with  the  persistence  of  symptoms  long  after  the  initiating  event.  This  same 
characteristic  may  also  explain  why  these  illnesses  appear  in  many  ways  resistant  to  treatment. 

4.  An  early  model  of  neuroinflammation  (extension  to  Task  3).  Elevated  levels  of  pro-inflammatory 
cytokines  negatively  impact  learning,  memory  and  neurogenesis.  The  intense  immune  activation  in  the  brain 
that  characterizes  infections,  injury,  neurotrauma  and  severe/chronic  stressful  conditions,  can  induce  hyper¬ 
excitability  of  neuronal  circuits  perpetuating  an  inflammatory  state  within  the  CNS  resulting  in  excitotoxicity, 
and  eventually  apoptosis  and  neurodegeneration  resulting  in  learning  and  memory  impairments.  To 
explore  these  mechanisms  we  have  constructed  a  first  model  with  Neurons,  Neural  Progenitor  Cells 
(NPCs),  Endothelial  Cells  (ECells),  Microglia,  and  Astrocytes  as  key  cellular  components,  while  interleukin 
(IL)-1,  IL-4,  IL-6,  tumor  necrosis  factor  (TNF)-a,  and  OX-2  membrane  glycoprotein  (CD200)  comprise  a 
simplified  neuro-inflammatory  response.  As  numerous  studies  show  communication  of  inflammatory 
information  to  the  brain  via  both  humoral  and  neuronal  mechanisms,  hormone  signaling  is  included  via 
Insulin-like  Growth  Factor  1  (IGF-1),  Vascular  Endothelial  Growth  Factor  (VEGF),  Brain  Derived 


Neurotrophic  Factor  (BDNF),  and  cortisol  (CORT).  The  neurotransmission  component  is  conveyed  with  the 
inclusion  of  Acetylcholine  (ACh),  Norepinephrine  (NE),  Glutamate  (Glut),  and  Adenosine  Triphosphate 
(ATP).  Stress-induced  immune  activation  is  a  neurally  initiated  phenomenon,  via  the  activation  of 
noradreneregic  pathways  and  altered  cholinergic  neurotransmission.  Elevated  brain  cytokines  produce 
further  activation  of  stress  response  systems  such  as  the  HPA  axis  and  the  SNS.  The  multiple  feed  forward 
and  feedback  connections  between  these  elements  found  in  the  neurophysiology  literature  are  depicted  in 
the  circuit  model  shown  in  Figure  5. 

•  Stable  immune  response  modes  in  the  brain  -  a  first  analysis:  Interaction  among  these 
various  cell  populations  via  immune,  hormone  and  neurotransmitter  signals  ultimately  revealed 
two  steady  states.  The  first,  again,  is  the  normal  reference  state  of  health  (SSO).  The  alternate 
steady  state  (SSI)  is  characterized  by  low  levels  of  ACh,  BDNF,  IGF-1 ,  IL-4,  and  VEGF  and 
suppressed  activity  of  Astrocytes,  ECells,  NPCs,  and  Neurons,  accompanied  by  elevated  levels 
of  CORT,  Glut,  IL-1 ,  IL-6,  and  TNF-a,  and  over  activation  of  Microglia.  This  is  consistent  with  a 
chronic  neuroinflammatory  state.  The  elevated  CORT  levels,  seen  to  align  with  GWI  in  our  other 
models,  suggests  a  possible  involvement  of  a  persistent  and  stable  neuro-inflammatory  cascade 
in  this  illness. 

5.  Continued  development  of  treatment  design  (Task  6,  7).  Analysis  of  the  above-mentioned  regulatory 
signaling  circuits  not  only  provides  information  describing  the  stable  steady  states  available  to  the  system 
but  also  extensively  describes  the  ensemble  of  transitory  states  that  lead  unequivocally  to  one  steady  state 
or  another;  these  are  said  to  lie  within  that  steady  state’s  basin  of  attraction.  Importantly,  these  subsets  of 
transitory  states  will  lead  to  that  specific  stable  state  independently  of  an  individual’s  immune  and  endocrine 
response  kinetics.  This  guaranteed  convergence  to  a  healthy  equilibrium  makes  them  attractive  as  broadly 
applicable  treatment  destination  states.  In  the  design  of  minimally  invasive  interventions  our  basic  paradigm 
is  therefore  to  identify  the  closest  transitory  state(s)  that  lie  within  the  basin  of  attraction  that  ensures  a 
return  to  normal  homeostasis. 

o  A  global  trajectory  search  formalism.  We  have  formalized  the  treatment  course  as  a  vector 
describing  a  path  from  the  state  of  disease  to  health.  Allowable  transitions  between  states  along 
the  path  consist  of  normal  evolution  of  the  system,  as  described  by  our  logic  rules,  and 
transitions  induced  by  clinically  feasible  interventions.  To  find  treatment  course  paths  that  meet 
these  criteria  we  have  formulated  our  search  as  a  global  optimization  problem.  We  have  chosen 
to  use  a  Genetic  Algorithm  (GA)  optimization  method,  as  the  discrete  nature  of  our  model 
naturally  accommodates  the  GA  solution  procedure.  Initially,  the  GA  seeds  the  solution  space  by 
generating  random  solutions  composed  of  binary  strings  or  “chromosomes”  representing  a 
treatment  path.  Each  member  of  this  initial  population  is  checked  against  a  fitness  function  and 
assigned  a  fitness  score.  Top  ranking  members  of  the  population  are  then  chosen  as  the  parent 
solutions  for  the  next  generation.  Each  generation  is  made  up  of  the  chosen  parent  population 
plus  combinations  of  crossed-over  “mated”  parent  solutions  with  a  small  chance  for  random 
mutation.  This  process  runs  over  a  set  number  of  generations  or  until  optimum  results  are  found. 
This  allows  a  rapid  search  of  the  global  space,  while  mutations  minimize  the  chance  of  remaining 
in  local  minima. 

o  A  multiple  objective  criterion.  Our  fitness  function  divides  the  overall  the  solution  string  into 
segments  describing  each  time-step  in  a  treatment  course.  The  overall  desirability  of  a  solution  is 
assessed  on  the  basis  of  three  objectives:  (i)  feasibility  or  compliance  with  the  model,  (ii) 
compliance  with  allowable  treatment  perturbations,  and  (iii)  the  minimal  invasiveness  and 
duration  of  treatment.  The  first  and  second  of  these  objectives  has  been  implemented  in  the  first 
trial  version.  For  each  time  step  segment  subsequent  time-steps  are  compared  to  the  allowable 
transition  states  and  assigned  a  compliance  score  based  on  the  minimum  hamming  distance 
separating  the  proposed  solution  state  and  the  allowable  states.  Overall  fitness  of  a  solution  is 
then  the  sum  of  the  hamming  measures  for  all  state  transitions  along  a  given  solution  path.  A 
fitness  value  of  zero  indicates  a  perfect  compliance  with  model  behavior  and  allowable 
interventions. 
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Based  on  this  paradigm  we  have  started  work  on  Tasks  6  and  7.  A  first  fitness  function  based  on  model 
compliance,  as  well  as  the  GA  algorithm  itself,  have  been  designed  and  implemented.  Currently,  we  are  in 
the  process  of  optimizing  code  parameters  (population  size,  number  of  generations,  mutation  rate,  etc...)  to 
efficiently  search  the  large  state  spaces  of  our  multi-systems  model. 

Additionally,  we  are  currently  refining  these  models  as  well  as  the  treatment  search  algorithm  to  incorporate 
the  effects  of  timescale.  This  will  make  it  possible  to  take  advantage  of  saddle  point  states  or  unstable 
intermediate  states  that  lie  between  the  basins  of  attraction.  We  are  currently  investigating  avenues  for 
exploiting  broad  classes  of  kinetic  scales  that  might  make  it  possible  to  reduce  the  treatment  complexity 
even  further  and  tailor  these  interventions  to  patient  sub-groups. 

6.  Continuing  work.  Ongoing  work  involves  the  continued  refinement  of  a  circuit  model  describing 
mechanisms  of  neuroinflammation  and  neurotransmission  in  the  brain.  Efforts  are  also  now  shifting  to  the 
completion  and  validation  of  the  treatment  design  algorithm.  This  will  be  the  major  area  of  development  as 
we  begin  simulation  of  treatment  strategies,  the  principal  deliverable  of  this  project. 

Timeline.  As  described  in  the  previous  report  dated  September  30,  2011,  the  University  of  Alberta’s 
Research  services  Office  submitted  on  behalf  of  the  principal  investigator  a  request  for  a  one-year 
extension  of  the  project  term  due  to  administrative  delays.  This  request  was  reviewed  initially  by  Ms.  Strock 
and  Dr.  Phillips  of  the  DoD  (January  23,  201 2)  and  we  were  asked  to  resubmit  this  request  at  a  later  date 
(6-8  months  before  end  of  project  term).  We  have  since  confirmed  with  Dr.  Rebecca  Fisher  that  this 
continues  to  be  the  correct  course  of  action  (ref.  email  from  Dr.  Fisher  dated  September  21, 2012).  In 
accordance  with  Dr.  Fisher’s  recommendation  we  are  submitting  a  formal  request  for  a  one-year  no-cost 
extension  as  part  of  our  request  to  transfer  this  award  to  Nova  Southeastern  University  retroactive  to  June 
1,  2013. 

Key  Research  Accomplishments. 

In  keeping  with  the  milestones  described  in  the  project  submission  initial  efforts  were  directed  at: 

•  Consistent  with  the  previously  completed  Task  1 ,  we  have  further  improved  the  efficiency  of  the 
serial  C  code,  again  delivering  order  of  magnitude  improvements  in  execution  speed  and  memory 
usage.  Importantly  we  have  engineered  a  parallel  framework  based  on  MPI  and  Pthread  protocols  to 
deploy  this  code  onto  distributed  high-performance  platforms.  This  code  is  now  deployed  and  fully 
operational  on  the  University  of  Miami  Pegasus  2  platform. 

•  Consistent  with  the  now  completed  Task  2,  we  have  continued  to  refine  our  previous  model  of 
immune  signaling  mechanisms.  These  now  include  the  actions  of  Th3  and  Th17  axes,  implemented 
in  models  at  two  levels  of  granularity. 

•  We  have  now  basically  completed  Task  3  as  defined  originally.  In  this  regard  we  have  produced  a 
refined  multi-axis  model,  further  developed  our  validation  scheme  and  submitted  to  PLoS  One  a  first 
complete  manuscript  describing  co-regulation  across  HPA,  HPG  and  immune  axes  in  men  and 
women. 

•  In  an  extension  to  original  Task  2  and  3,  we  have  produced  a  first  circuit  model  of  inflammatory 
processes  occurring  in  the  brain  and  involving  the  cell  types  and  immune  signaling  specific  to  this 
physiological  compartment.  Early  analyses  of  this  model  show  the  persistence  of  a  chronic 
neuroinflammatory  state  perpetuated  by  overactive  microglia  and  underactive  astrocytes,  leading  to 
loss  of  neuron  function,  in  conjunction  with  elevated  levels  of  cortisol. 

•  Consistent  with  Task  4  and  Task  5  we  have  conducted  a  refined  analysis  of  multi-stability  properties 
of  both  the  broad  HPA-HPG-immune  model  and  the  detailed  BIS  immune  model.  Comparing 
predicted  equilibrium  states  with  experimental  immune  and  endocrine  data  from  male  and  female 
GWI  and  CFS  subjects  we  find: 

o  Male  GWI  and  CFS  subjects  align  with  states  showing  hypercortisolism,  low  testosterone, 
elevated  Thl  inflammatory  cytokines,  decreased  NK  cell  activity  in  conjunction  with  an 
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elevated  Th17  response,  although  male  CFS  subjects  also  show  a  propensity  to  align  with  an 
elevated  Treg  response  suggesting  a  mixed  immune  signature  for  this  illness  not  seen  in 
GWI. 

o  Female  CFS  subjects  align  with  states  showing  hypocortisolism,  elevated  estrogen,  and  a 
shift  towards  Th2  activation. 

•  We  have  re-assessed  our  approach  to  treatment  design  (Task  6,  7)  and  have  begun  encoding  an 
approach  based  on  a  global  search  for  a  treatment  course  assisting  an  optimal  walk  through  a 
discrete  endocrine  immune  state  space  leading  from  an  illness  to  a  healthy  condition. 

Reportable  Outcomes. 

The  results  of  these  latest  analyses  are  being  communicated  as  follows: 

•  The  previous  draft  manuscript  Craddock  et  al.,  2013,  enclosed  as  Appendix  A,  has  been  extensively 
revised  and  is  now  submitted  to  the  journal  PLoS  One.  Similarly  we  expect  the  extensions  and 
revisions  to  the  detailed  immune  model  (working  document  in  Appendix  B,  Annual  Report  2012, 
Fritsch  et  al.,  2013),  to  be  ready  for  submission  to  the  journal  Molecular  Systems  Biology  by  October 
this  year. 

•  Early  results  were  presented  at  a  closed  meeting  sponsored  by  the  CDC  and  the  CFIDS  Association 
of  America  and  held  at  the  Cold  Spring  Harbor  Laboratory’s  Banbury  Centre  in  Long  Island,  NY 
(Sep.  30 -Oct  3,  2012). 

•  We  will  be  submitting  two  abstracts  for  oral  presentation  at  the  IACFS/ME  1 1th  Biennial  International 
Research  and  Clinical  Conference  to  be  held  in  San  Francisco,  California,  USA,  March  20-23,  2014. 
The  conference  is  co-sponsored  by  Stanford  University. 

Regarding  synergy  with  complementary  research  efforts,  these  findings  were  recently  used  to  secure  an 
invited  GWIRP  Consortium  Award,  now  awarded  (prime  institution  -  Nova  Southeastern  University).  The  are 
also  being  used  in  support  of  2  VA  Merit  applications  that  have  been  reviewed  and  invited  for  resubmission 
this  month. 

Conclusions. 

We  are  currently  processing  a  formal  request  for  a  one-year  no-cost  extension  of  the  project  term  due  to  a 
delayed  start.  We  have  carried  out  a  major  shift  in  paradigm  and  continue  to  refine  these  regulatory  circuit 
models  as  well  as  developing  new  components  such  as  the  neuroinflammatory  model.  The  basic 
algorithmic  framework  has  now  been  translated  and  re-engineered  to  deploy  larger  more  detailed  models 
on  distributed  high-performance  platforms  like  the  Pegasus  2  platform  at  the  University  of  Miami. 

Simulations  based  on  these  models  have  shown  that  the  illness-specific  effects  of  gender  are  particularly 
striking.  Work  continues  on  the  refinement  of  the  intervention  design  component.  Initial  analyses  favor  the 
deployment  of  a  joint  hormone-immune  intervention  over  strategies  that  target  these  systems  separately. 

Personnel  receiving  support  from  this  award: 

•  Travis  Craddock,  Ph.D.,  Senior  Research  Associate,  Broderick  Clinical  Systems  Biology  Laboratory, 
University  of  Alberta; 

Now  co-investigator  and  Assistant  Professor  of  Psychology,  Nova  Southeastern  University,  FL. 
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•  Mark  Rice,  undergraduate  student  computer  science,  Research  Intern  and  Chief  Programmer,  Nova 
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Version  Performace  Measures:  Run  Time  vs.  Ternary  Nodes 


Ternary  Nodes 


—Version  0 

- Version  1 

- Version  2 

Version  3m0 

- Version  3ml 

Version  3m2 
—Version  3m3 
- Version  3m4 


Figure  1.  Evolution  of  computational  performance.  Evolution  of  computer  wall  time  as  a  function  model 
complexity  described  in  terms  of  the  number  of  state  variables  (ternary  nodes).  In  less  than  a  year,  re¬ 
engineering  of  the  computer  code  supporting  the  identification  of  stable  states  in  a  regulatory  system  has 
enable  an  almost  2-fold  increase  in  the  number  of  state  variables  in  the  circuit  model. 
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F  G 


Figure  2.  Increased  granularity  multi-axis  model.  A  significant  revision  of  the  discrete  circuit  model  of  HPA 
function  (A,B)  augmented  with  HPG-immune  interactions  in  male  subjects  (C)  and  female  subjects  (D-G)  in  the 
specific  case  of  positive  feedback  along  the  female  HPG  axis  and  suppressive  interaction  with  the  HPA  axis 
(revised  and  resubmitted  manuscript)  [6].  Green  directed  edges  represent  an  up-regulation  of  the  target  by  the 
source  node  whereas  a  red  terminal  edge  represents  a  suppressive  action. 
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Figure  3.  Detailed  Immune  model  revised.  Circuit  diagram  of  the  detailed  immune  system  model  revised  to 
include  elements  of  Th17  and  Treg  activity  mediated  by  TGF-p,  IL-21,  IL-23,  IL-27  and  others.  This  is  a 
significant  increase  in  granularity  from  the  previous  such  model  and  has  resulted  in  a  revision  of  draft 
manuscript,  now  underway  and  due  for  submission  before  year  end  [12], 
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Figure  4.  High-resolution  immune  model  with  neurotransmission  oversight.  Circuit  diagram  of  a  first  prototype 
model  capturing  fine-grained  immune  signaling  with  the  contribution  of  immune  modulating  neurotransmitters 
neuro-peptide  Y  (NPY),  acetylcholine  (ACh)  and  norepinephrine  (NEpi).  This  model  is  still  in  progress  and  will 
incorporate  the  Th17  and  Treg  axes  as  well  as  additional  neurotransmitters  as  we  move  forward. 
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Figure  5.  A  first  prototype  model  of  brain  immunity.  A  first  circuit  model  describing  the  regulation  of  neuro¬ 
inflammation  in  the  brain  that  includes  the  role  of  neurons,  neural  progenitor  cells  (NPCs),  endothelial  cells 
(ECells),  microglia,  and  astrocytes  as  key  cellular  components  as  well  as  basic  signaling  mechanisms  involving 
interleukin  (IL)-1,  IL-4,  IL-6,  tumor  necrosis  factor  (TNF)-a,  and  OX-2  membrane  glycoprotein  (CD200)  and 
other  molecular  messengers. 
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Appendix  A: 

Craddock  TJA,  Fritsch  P,  Rice  MA  Jr,  del  Rosario  RM,  Miller  DB,  Fletcher  MA,  Klimas  NG,  Broderick  G.  A  Role 
for  Homeostatic  Drive  in  the  Perpetuation  of  Complex  Chronic  Illness:  Gulf  War  Illness  and  Chronic  Fatigue 
Syndrome.  2013,  Submitted  PLoS  One,  Under  editorial  review. 
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27  Abstract 

28  A  key  component  in  the  body’s  stress  response,  the  hypothalamic-pituitary-adrenal  (HPA) 

29  axis  orchestrates  changes  across  a  broad  range  of  major  biological  systems.  Its  dysfunction 

30  has  been  associated  with  numerous  chronic  diseases  including  Gulf  War  Illness  (GWI)  and 

3 1  chronic  fatigue  syndrome  (CFS).  Though  tightly  coupled  with  other  components  of  endocrine 

32  and  immune  function,  few  models  of  HPA  function  account  for  these  interactions.  Here  we 

33  extend  conventional  models  of  HPA  function  by  including  feed-forward  and  feedback 

34  interaction  with  sex  hormone  regulation  and  immune  response.  We  use  this  multi-axis  model 

35  to  explore  the  role  of  homeostatic  regulation  in  perpetuating  chronic  conditions,  specifically 

36  GWI  and  CFS.  An  important  obstacle  in  building  these  models  remains  the  scarcity  of  in  vivo 

37  kinetic  data.  We  circumvented  this  using  a  discrete  logic  representation  based  solely  on 

38  literature  of  physiological  and  biochemical  connectivity  to  provide  a  qualitative  description  of 

39  system  behavior.  This  connectivity  model  linked  molecular  variables  across  the  HPA  axis, 

40  hypothalamic-pituitary-gonadal  (HPG)  axis  in  men  and  women,  as  well  as  a  simple  immune 

41  network.  Inclusion  of  these  interactions  produced  at  multiple  alternate  homeostatic  states. 

42  Experimental  data  for  endocrine-immune  markers  measured  in  male  GWI  subjects  showed 

43  the  greatest  alignment  with  predictions  of  a  naturally  occurring  alternate  steady  state 

44  presenting  with  hypercortisolism,  low  testosterone  and  a  shift  towards  a  Thl  immune 

45  response.  In  female  CFS  subjects,  expression  of  these  markers  aligned  with  an  alternate 

46  homeostatic  state  displaying  hypocortisolism,  high  estradiol,  and  a  shift  towards  an  anti- 

47  inflammatory  Th2  activation.  These  results  support  a  role  for  homeostatic  drive  in 

48  perpetuating  dysfunctional  cortisol  levels  through  persistent  interaction  with  the  immune 

49  system  and  HPG  axis.  This  same  basic  drive  may  also  perpetuate  sexually  dimorphic 

50  responses  due  to  inherently  different  behavior  of  the  male  and  female  HPG.  Though  coarse, 

51  these  models  may  nonetheless  support  the  design  of  robust  treatments  that  might  exploit 

52  these  regulatory  regimes. 

53 
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53  Introduction 

54  The  hypothalamic-pituitary-adrenal  (HPA)  axis,  a  key  component  in  the  body’s  stress 

55  response,  serves  to  articulate  changes  in  a  broad  range  of  homeostatic  regulators  as  a 

56  function  of  environmental  cues.  Such  cues  can  consist  of  both  physical  stressors  (injury, 

57  infection,  thermal  exposure)  and  psycho-emotional  stressors  (frustration,  fear,  fight  or  flight 

58  decisions).  Instantiation  of  this  survival  program  is  accomplished  through  controlled 

59  modulation  of  the  neuroendocrine  and  immune  systems,  as  well  as  the  sympathetic  nervous 

60  systems  [1-3].  Considering  its  function  as  a  broad-reaching  integrator  of  major  physiological 

61  systems,  it  is  no  surprise  that  numerous  chronic  conditions  have  been  associated  with 

62  abnormal  regulation  of  the  HPA  axis,  including  major  depressive  disorder  (MDD)  [4,  5],  post- 

63  traumatic  stress  disorder  (PTSD)  [6-8],  Alzheimer’s  disease  [9],  Gulf  War  Illness  (GWI)  [IQ- 

64  12],  and  chronic  fatigue  syndrome  (CFS)  [13-15].  When  compared  to  non-deployed 

65  veterans,  Golier  et  al.  [1 0]  found  that  symptomatic  Gulf  War  veterans  without  psychiatric 

66  illness,  as  well  as  veterans  with  PTSD  alone,  showed  significantly  greater  cortisol 

67  suppression  to  dexamethasone  (DEX)  suggesting  markedly  enhanced  negative  feedback 

68  along  the  HPA  axis.  Further  study  by  these  same  investigators  indicated  that  this  might  be 

69  due  to  a  significantly  attenuated  ACTH  response  by  the  pituitary  in  veterans  with  GWI 

70  without  PTSD  [11,  12],  A  similar  suppression  of  cortisol  response  to  DEX  was  found  in  CFS 

71  subjects  by  Van  Den  Eede  et  al.  [13]  with  this  being  further  exacerbated  by  oestrogen 

72  intake.  With  regard  to  HPA  circadian  dynamics,  CFS  subjects  were  found  to  exhibit 

73  significantly  increased  adrenal  sensitivity  to  ACTH  and  marginally  increased  inhibitory 

74  feedback  during  the  nocturnal  period  when  compared  with  control  subjects  and  CFS 

75  subjects  comorbid  with  fibromyalgia  (FM)  [14,  15],  Conversely  the  pain-dominant  CFS-FM 

76  subjects  showed  significantly  blunted  cortisol  inhibitory  feedback.  While  evidence  such  as 

77  this  implicates  abnormal  regulation  of  HPA  function  leading  to  chronic  hypocortisolic  and 

78  hypercortisolic  states  in  these  illnesses,  the  genesis  of  this  dysregulation  is  unclear. 

79 


80  Previously  we  investigated  the  possibility  that  some  of  these  pathological  states  may 

81  coincide  with  naturally  occurring  alternate  homeostatic  stable  states  [16].  These  “backup 

82  programs”  would  offer  a  way  of  maintaining  homeostatic  control  in  crisis  situations  at  the 

83  cost  of  reduced  function.  The  existence  of  such  multiple  stable  states  is  characteristic  of 

84  systems  that  incorporate  feed-forward  and  feedback  mechanisms.  Feedforward  loops  in 

85  biology  play  the  crucial  role  of  driving  rapid  acute  responses,  while  feedback  loops  will 

86  generally  limit  the  extent  of  a  response.  Both  will  also  drive  complex  dynamic  behavior, 

87  including  differentiation  and  periodicity  [17],  While  small  perturbations  may  force  temporary 

88  departures,  these  systems  return  to  their  original  resting  states  once  these  perturbations  are 

89  removed.  If  however,  the  perturbation  is  of  significant  strength  and  duration,  the  system 

90  may  be  incapable  of  returning  to  its  normal  operating  regime  and  instead  may  assume  a 

91  new  alternate  resting  state.  Knowledge  of  the  system  dynamics  can  allow  us  to  map  these 

92  different  stable  states  and  several  mathematical  models  of  the  HPA  exist  [18-26],  So  far, 

93  only  one  such  model  is  known  to  accommodate  multi-stability  in  the  dynamic  behavior  of  the 

94  HPA  axis.  It  does  so  via  the  addition  of  a  feed-forward  mechanism  involving  dimerization  of 

95  the  glucocorticoid  receptor  (GR)  complex  [27]  (Figure  1).  In  this  process  glucocorticoid 

96  (GC)  bound  GRs  form  homodimers  that  translocate  into  the  cell  nucleus  to  bind  DNA,  up- 

97  regulating  GR  synthesis  and  producing  a  positive  feedback  loop.  However,  this  model  and 

98  the  majority  of  other  models  do  not  extend  beyond  the  physiological  boundaries  of  the  HPA 

99  axis  itself  and  thus  are  limited  in  their  predictive  capabilities.  As  discussed  in  the  following 

100  sections,  HPA  activity  is  intertwined  with  the  behavior  of  the  hypothalamic-pituitary-gonadal 

101  (HPG)  axis  and  the  immune  system,  among  others,  and  this  interplay  should  not  be  ignored 

102  when  considering  the  number  and  nature  of  stationary  states  available  to  the  overarching 

103  system.  Our  hypothesis  is  that  these  alternate  regulatory  regimes  may  facilitate  the 

104  persistence  of  complex  chronic  illnesses  like  GWI  and  CFS.  To  evaluate  the  role  of  alternate 

105  homeostatic  attractors  in  these  illnesses  we  constructed  a  computational  model  of  regulatory 

106  control  linking  the  HPA,  HPG  and  immune  systems. 

107 


A 


108  There  is  a  substantial  body  of  physiological  and  biochemical  data  for  many  biological 

109  systems  describing  the  connectivity  between  molecular  and  cellular  elements,  the  presence 

1 10  of  recurring  structural  motifs  and  functional  modules.  For  example,  negative  autoregulation, 

111  in  which  a  transcription  factor  represses  its  own  transcription,  is  a  simple  network  motif 

112  observed  in  many  transcription  networks.  While,  numerous  motifs  have  been  found  in 

1 13  biological  networks  (negative/positive  autoregulation,  coherent/incoherent  and  multi-output 

114  feedforward  loops,  single-input  modules  and  dense  overlapping  regulons)  [28],  data 

115  regarding  the  precise  stoichiometry  and  kinetics  of  these  systems  in  humans  is  extremely 

1 16  limited.  Many  existing  models  rely  heavily  on  animal  data  as  a  source  of  kinetic  parameters, 

1 17  or  adopt  general  order  of  magnitude  estimates  when  this  data  is  lacking.  To  circumvent  this 

118  issue  and  draw  on  the  rich  body  of  known  molelcular  and  cellular  interactions  in 

1 19  physiological  and  biochemistry,  we  have  adopted  the  discrete  logical  network  methodology 

120  proposed  originally  by  Thomas  et  al.  [29,  30]  and  developed  further  by  Mendoza  and 

121  Xenarios  [31],  By  applying  logic  rules  to  a  network  of  known  interactions  it  is  possible  to 

122  identify  the  number  of  stable  resting  states,  their  type  as  well  as  their  molecular  and  cellular 

123  description,  without  detailed  knowledge  of  the  response  dynamics.  In  this  work  we  use  this 

124  method  to  extend  our  previous  analysis  of  human  HPA  axis  dynamics  by  including  its 

125  regulatory  interactions  with  the  neighboring  HPG  axis  and  immune  system.  This  resulting 

126  mathematical  model  better  represents  the  complexity  of  endocrine-immune  interactions  by 

127  supporting  the  detection  and  identification  of  alternate  resting  modes  of  the  HPA-HPG- 

128  immune  axis.  Based  on  connectivity  information  alone,  we  show  that  multi-stability  is  easily 

129  obtained  from  these  interacting  systems.  Moreover,  we  show  that  experimental  data  from 

130  our  on-going  studies  of  GWI  and  CFS  show  better  alignment  with  these  alternate  resting 

131  modes  than  with  the  typical  healthy  homeostatic  stable  state.  Ultimately,  knowledge  of  such 

132  homeostatic  modes  could  be  used  to  identify  promising  applications  of  pharmaceutical, 

133  hormone  and/or  immune  therapy  that  exploit  the  body’s  natural  dynamics  to  reinforce 

134  treatment  effects. 

135 


136  Methods 

137  Ethics  Statement 

138  All  subjects  signed  an  informed  consent  approved  by  the  Institutional  Review  Board  of  the 

139  University  of  Miami.  Ethics  review  and  approval  for  data  analysis  was  also  obtained  by  the 

140  IRB  of  the  University  of  Alberta. 

141 

142  An  Integrative  Multi-systems  Model  of  the  HPA-HPG-lmmune  System 

143  There  is  a  substantial  amount  of  physiological  data  describing  the  HPA,  HPG  and  immune 

144  systems  as  stand-alone  entities.  To  a  much  lesser  degree  there  also  exists  evidence  for  the 

145  mutual  interactions  between  these  systems.  The  following  sections  describe  the 

146  experimental  evidence  used  to  infer  the  topology  of  an  overarching  HPA-HPG-immune 

147  interaction  network  (Figure  1). 

148 

149  The  HPA  Axis:  Activation  of  the  HPA  axis  begins  at  the  paraventricular  nucleus  (PVN)  of 

150  the  hypothalamus.  Specifically,  afferents  transmitting  stress  related  information  in  the  brain 

151  converge  on  the  medial  parvocellular  neurons  of  the  PVN  inducing  the  release  of  several 

152  peptides,  including  corticotropin-releasing  hormone  (CRH)  and  arginine  vasopression  (AVP), 

153  into  the  pituitary  hypophysial-portal  circulation.  The  unique  vascular  system  allows  very 

154  small  quantities  of  these  hypothalamic  hormones  to  act  directly  on  their  targets  in  the 

155  anterior  pituitary  without  dilution  by  systemic  circulation.  CRH  and  AVP  act  in  conjunction  on 

156  membrane  bound  CRH-R1  receptors  in  the  anterior  pituitary  to  stimulate  adrenocorticotropic 

157  hormone  (ACTH)  synthesis,  and  its  rapid  release  into  peripheral  circulation.  ACTH 

158  circulates  to  the  adrenal  cortex  where  it  acts  on  the  membrane  bound  MC2-R  receptor  to 

159  simulate  the  release  of  GCs  (corticosterone  in  the  rat,  and  cortisol  (CORT)  in  humans  and 

160  nonhuman  primates).  To  regulate  the  stress  response,  GCs  exert  negative  feedback  at  the 

161  hypothalamus  and  pituitary  to  inhibit  further  synthesis  and  release  of  CRH  and  ACTH, 

162  respectively  [32],  This  is  the  standard  view  of  the  HPA  axis  utilized  in  the  majority  of  models 

163  (Figure  1  A).  However,  as  noted  by  Gupta  et  al.  [27]  circulating  glucocorticoids  act  via 
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164  cytostolic  GRs,  which,  unlike  membrane  bound  receptors,  dimerize  (GRD)  and  translocate 

165  into  the  cell  nucleus  upon  activation  to  up-regulate  GR  synthesis  and  interact  with  other 

166  relevant  transcription  factors,  or  GC-sensitive  genes  (Figure  1  B).  Gupta  et  al.  included  this 

167  GR  expression  feedforward  loop  at  the  pituitary,  as  it  is  a  main  driver  of  the  HPA  axis,  and 

168  found  a  resulting  bistability  in  the  HPA  system  [21].  However,  all  nucleated  cells  possess 

169  GRs,  as  GCs  influence  practically  every  system  in  the  body,  suggesting  this  feedforward 

170  loop  may  be  important  in  other  tissues  beyond  the  HPA  axis.  As  described  below  major 

171  systems  affected  by  GCs  include  the  HPG  axis  and  immune  system. 

172 

173  The  HPG  Axis:  GCs  have  an  inhibitory  effect  on  the  HPG  axis,  a  central  regulator  of  the 

174  reproductive  system,  at  all  levels  [33-37],  Activation  of  the  HPG  starts  from  brain  generated 

175  pulsatile  signals  that  stimulate  the  preoptic  area  of  the  hypothalamus  to  produce 

176  gonadotropin-releasing  hormone  (GnRH).  GnRH  is  secreted  into  the  pituitary  hypophysial 

177  portal  bloodstream,  which  carries  it  to  the  pituitary  gland,  where  it  activates  membrane 

178  bound  GnRH-R  receptors,  resulting  in  the  synthesis  and  secretion  of  luteinizing  homone 

179  (LH)  and  follicle-stimulating  hormone  (FSH)  into  circulation.  These  gonadotropins  flow  to  the 

180  gonads  where  they  work  synergistically  to  promote  the  secretion  of  the  sex  steroids.  In 

181  males,  LH  binds  to  receptors  on  Leydig  cells  in  the  testes  to  stimulate  the  synthesis  and 

182  secretion  of  testosterone  (TEST).  In  females,  LH  activates  receptors  on  Theca  interna  cells 

183  in  the  ovaries  to  stimulate  the  release  of  androstenedione,  which  is  aromatized  by  granulosa 

184  cells  to  produce  estradiol  (EST),  and  progesterone  (PROG).  TEST  negatively  feeds  back  on 

185  the  HPG  to  inhibit  GnRH,  FSH  and  LH  secretion  and  synthesis  [33],  This  feedback 

186  mechanism  is  somewhat  more  complex  in  females  where,  depending  on  the  phase  of  the 

187  female  menstrual  cycle,  EST  and  PROG  can  exert  either  positive  or  negative  feedback  on 

188  the  production  and  release  of  GnRH  and  the  gonadotropins  [36,  38,  39], 

189 

190  A  lesser-known  aspect  is  that  several  components  of  the  HPG  axis  exert  reciprocal  effects 

191  on  the  HPA  axis  [33,  34,  36].  Testosterone  exhibits  an  inhibitory  effect  on  all  levels  of  the 
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192  HPA  [33]  (Figure  1  C),  whereas  EST  and  PROG  can  serve  to  stimulate  or  inhibit  the  HPA 

193  axis  depending  on  menstrual  cycle  phase,  or  phase  of  life  [34].  These  affects  may  be 

194  mediated  through  changes  in  adrenocorticoid  synthesis,  stress-induced  ACTH  and  GC 

195  release,  and  CRH  and  AVP  synthesis  in  the  PVN,  by  direct  activation  of  oestrogen  and 

196  androgen  receptors  along  the  HPA  or  via  interaction  between  GRs  and  sex  steroid  receptors 

197  to  regulate  transcription  [33,34,36],  Thus,  an  interactive  functional  crosstalk  exists  between 

198  the  HPA  and  HPG  axes,  which  cannot  be  ignored  when  investigating  HPA  axis  regulation 

199  and  dysfunction.  Mutual  inhibition  between  the  HPA  and  HPG  (Figure  1  C)  was  considered 

200  standard  for  males.  However,  as  it  is  not  clear  whether  the  EST  and  PROG 

201  inhibition/stimulation  of  the  HPA  occurs  in  coordination  with  the  inhibition/stimulation  of  the 

202  HPG,  these  cases  were  explored  for  females  alone  as  separate  alternative  models  of  the 

203  HPA-HPG  interaction  (Figure  1  D-G)  in  addition  to  the  model  considered  for  males. 

204 

205  A  Simple  Model  of  the  Immune  System:  While  not  typically  considered  part  of  the 

206  neuroendocrine  system,  the  immune  system  plays  a  very  important  role  in  regulating  the 

207  HPA  axis.  Here  we  base  our  simplified  immune  system  upon  our  previous  work  detailing  the 

208  communication  network  of  the  immune  response  [40],  Cells  of  the  innate  immune  response 

209  (ICells),  including  mononuclear  phagocytes,  such  as  macrophages,  and  dendritic  cells, 

210  natural-killer  (NK)  cells,  endothelial  cells  and  mucosal  epithelial  cells,  communicate  via  the 

211  release  of  numerous  cytokines.  Cytokines  that  regulate  the  innate  immune  response  (MR) 

212  include  interleukin  (IL)  -1 ,  IL-6,  IL-8  and  tumor  necrosis  factor  alpha  (TNF-a),  and  can  also 

213  include  IL-12,  a  primary  mediator  of  early  innate  immunity.  Primarily,  these  signals  serve  to 

214  activate  and  recruit  other  ICells,  which  in  turn  produce  more  cytokines.  IL-15,  which 

215  stimulates  proliferation  of  NK  cells  and  effector  T-lymphocytes,  can  also  be  considered  as 

216  part  of  the  MR  as  well  as  IL-23,  an  important  inflammatory  signal  contributing  to  the  Th17 

217  response  against  infection. 

218 
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219  MR  signals  can  also  serve  to  prime  helper  T  cells  towards  a  Thl  type  adaptive  immune 

220  response  (T 1  Cell).  This  response  produces  Thl  proinflammatory  cytokines  (TICyt) 

221  including  IL-2,  interferon-gamma  (IFN-y),  and  tumor  necrosis  factor  beta  (TNF-|3),  which 

222  further  activates  ICells,  while  suppressing  the  Th2  adaptive  immune  response  (T2Cell). 

223  The  T2Cell  is  responsible  for  the  production  of  the  Th2  anti-inflammatory  cytokines  (T2Cyt) 

224  IL-4,  IL-5,  IL-10  and  IL-13,  which  have  important  anti-inflammatory  and  immunosuppressive 

225  activities,  and  serve  to  inhibit  the  activity  of  T 1  Cell  and  ICells. 

226 

227  Cytokines  can  also  serve  as  mediators  between  the  immune  and  endocrine  systems. 

228  Between  the  HPA  and  the  immune  network  there  exists  a  mutual  crosstalk  [41-43]  (Figure  1 

229  C-G).  The  MR  and  T 1  Cell  cytokines  selected  here  serve  to  stimulate  the  HPA  axis  at  all 

230  levels  [41-43],  CORT,  in  turn,  acts  to  suppress  the  activity  of  ICells  (specifically  NK  cells 

23 1  [44],  and  DC  cells  [45]),  and  the  T 1  Cell  [46]  causing  a  shift  from  the  inflammatory  to  the  anti- 

232  inflammatory  response  [41, 42,  47],  The  interaction  between  the  HPG  and  the  immune 

233  system  is  complex  and  sexually  dimorphic,  and  is  still  an  active  field  of  research.  However, 

234  at  a  general  coarse  level  of  description  TEST  serves  to  stimulate  the  development  of  the 

235  Thl  response  [48]  (Figure  1  C),  whereas  EST  inhibits  the  Thl  response  causing  a  shift 

236  towards  the  Th2  anti-inflammtory  response  [48,49],  The  reciprocal  crosstalk  from  the 

237  immune  system  to  the  HPG  is  equally  intricate.  In  broad  terms  this  conversation  is 

238  communicated  via  T 1  Cyt.  Receptors  for  TNF-a  and  IFN-y  are  expressed  in  testicular  Leydig 

239  cells  and  there  is  evidence  that  these  cytokines  can  directly  inhibit  testosterone  production 

240  [50].  TNFa  also  decreases  the  release  of  GnRH  in  the  hypothalamus  and  LH  in  the  pituitary 

241  gland  in  both  males  [50]  and  females  [51]  eventually  leading  to  a  decrease  in  sex  steroid 

242  levels.  As  such,  we  model  the  TICyt  as  inhibiting  GnRH  and  LH/FSH  in  both  male  and 

243  female  models. 

244 

245 

246 


Q 


247  A  Discrete  State  Representation 

248  Following  the  methods  of  Thomas  et  al.  [29,  30],  and  more  recently  Mendoza  and  Xenarios 

249  [31],  the  neuroendocrine-immune  system  was  represented  as  a  connectivity  model 

250  consisting  of  interconnected  molecular  and  cellular  variables  with  three  discrete  states:  -1 

251  (inhibited),  0  (nominal)  and  1  (activated).  According  to  this  type  of  model  the  current  state  of 

252  all  variables  in  a  system  is  described  by  a  state  vector  x(t) ,  such  that: 

253  xt=xlt,x2t,...,xNt  (1) 

254  where  xN(t)is  the  state  of  the  Nth  variable  of  the  system  at  time  t.  The  image  vector  x(t  + 1) 

255  describes  the  preferred  state  towards  which  the  system  evolves  in  the  next  time  increment. 

256  The  state  value  of  the  image  vector  for  the  ith  variable  is  determined  from  its  current  state 

257  and  a  set  of  balanced  ternary  logic  statements  based  on  the  current  value  of  variable  and 

258  the  mode  of  action  (i.e.  activate  or  inhibit)  of  the  neighboring  input  variables.  These  logic 

259  statements  are  expressed  as  follows  (Eq.  2): 

260  xit+l=(xil  AtVxi2  At...xi jAt)V(xil  Tt\lxi2 Tt...xik  Tt)(xil  AtVxi2  At...xi jAt)~’(xil  TtVxi2 Tt...xik  Tt) 

261  (2) 
262 

263  where  the  V,  v,  and  -  symbols  are  ternary  HIGH/LOW  PASS,  OR  and  NOT  operators,  x* 

264  is  the  state  of  the  ith  variable’s  fh  activator,  x'k  is  the  state  of  the  ith  variable’s  k th  inhibitor.  The 

265  ternary  operators  given  in  Equation  (2)  are  described  in  further  detail  in  Supplementary 

266  Tables  1-  3.  The  first  entry  in  Equation  (2)  is  used  when  the  variable  possesses  both 

267  activators  and  inhibitors,  the  middle  when  the  variable  has  only  activators  and  last  when  the 

268  activator  has  only  inhibitors. 

269 

270  Applying  Equation  (2)  to  each  variable  in  the  model  for  the  mth  state  of  the  system,  xm(t) , 

271  defines  the  image  vector  xm(t  + 1)  for  that  state.  With  xm(t  + 1)  defined,  the  system  maybe 

272  updated  asynchronously  (allowing  only  one  variable  to  change  at  a  time)  following  the 
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273  generalized  logical  analysis  of  Thomas  et  al.  [29,  30].  According  to  this  method  the  ith 

274  variable  of  the  mth  state  vector  xm(t)  is  moved  one  step  towards  its  preferred  image 

275  xm(t  +  l)  (e.g.  If  xm(t)  =  -1  and  xm(t  + 1)  =1,then  x^t)  is  set  to  0).  Thus,  for  each  current 

276  state  of  the  system  there  are  potentially  several  subsequent  states  towards  which  it  may 

277  asynchronously  evolve. 

278 

279  The  number  of  states,  and  the  values  they  can  be  assigned,  determine  the  total  number  of 

280  states  available  to  the  model  system.  With  the  ternary  logic  used  here,  a  model  of  N 

281  variables  possesses  3N  states.  As  a  result,  the  number  of  states  increases  rapidly  as  new 

282  variables  are  added.  By  analyzing  all  possible  states  of  the  system  a  temporal  sequence  of 

283  states  may  be  discerned.  To  interpret  the  results,  each  state  of  the  system  can  be 

284  represented  as  an  element  in  a  graph.  The  evolution  from  one  state  to  a  subsequent  state 

285  can  be  represented  as  a  directed  edge  between  the  two  states  in  this  graph.  Representation 

286  of  the  state  trajectories  in  this  fashion  makes  it  possible  to  draw  on  the  concepts  and  tools  of 

287  graph  theory  for  analysis  of  the  system  dynamics.  Steady  states  are  defined  as  those  states 

288  for  which  the  image  vector  is  the  same  as  the  current  state  vector;  in  other  words  the  state 

289  possesses  an  out  degree  of  0. 

290 

291 

292  Comparison  to  Model 

293  GWI  Cohort  Sample  Collection:  Similar  cytokine  profiles  and  endocrine  measures  were 

294  obtained  as  part  of  a  larger  ongoing  study  of  27  GWI  and  29  HC  subjects  recruited  from  the 

295  Miami  Veterans  Administration  Medical  Center.  Subjects  were  male  with  an  average  age  of 

296  43  years  and  BMI  of  28.  Inclusion  criteria  was  derived  from  Fukuda  et  al.  [52],  and  consisted 

297  in  identifying  veterans  deployed  to  the  theater  of  operations  between  August  8,  1990  and 

298  July  31,  1991,  with  one  or  more  symptoms  present  after  6  months  from  at  least  2  of  the 

299  following:  fatigue;  mood  and  cognitive  complaints;  and  musculoskeletal  complaints.  Subjects 


300  were  in  good  health  prior  to  1990,  and  had  no  current  exclusionary  diagnoses  [53].  Use  of 

301  the  Fukuda  definition  in  GWI  is  supported  by  Collins  et  al.  [54],  Control  subjects  consisted  of 

302  gulf  war  era  sedentary  veterans  and  were  matched  to  GWI  subjects  by  age,  body  mass 

303  index  (BMI)  and  ethnicity.  Additional  details  regarding  this  cohort  and  the  laboratory  assays 

304  performed  are  available  in  Broderick  et  al.  [55], 

305 

306  CFS  Cohort  Sample  Collection:  Levels  of  cortisol  (CORT)  and  estradiol  (EST)  measured 

307  in  peripheral  blood  were  obtained  from  the  Wichita  Clinical  dataset  [56]  for  a  group  of  39 

308  female  CFS  subjects  and  37  Healthy  controls  (HCs)  with  an  average  age  of  52  years  and  an 

309  average  body  mass  index  (BMI)  of  29.  Additional  details  of  this  cohort  and  the  laboratory 

310  assays  performed  may  be  found  in  work  previously  reported  by  our  group  [57,  58].  Multiplex 

311  cytokine  profiles  were  obtained  in  plasma  from  a  separate  but  demographically  comparable 

3 12  cohort  of  40  female  CFS  subjects  and  a  group  of  59  healthy  female  matched  control 

313  subjects  studied  by  our  group  at  the  University  of  Miami  [59],  Average  age  in  this  cohort  was 

314  53  years  with  an  average  BMI  of  26.  Profiling  of  cytokine  concentrations  was  performed  in 

315  morning  blood  plasma  samples  using  an  enzyme-linked  immuno-absorbent  assay  (ELISA)- 

316  based  assay.  Details  of  this  protocol  and  results  of  a  comparative  analysis  of  cytokine 

317  expression  patterns  are  available  in  Broderick  et  al.  [59].  In  both  studies  a  diagnosis  of  CFS 

318  was  made  using  the  International  Case  Definition  [53,60].  Exclusion  criteria  for  CFS  included 

319  all  of  those  listed  in  the  current  Centers  for  Disease  Control  (CDC)  CFS  case  definition,  as 

320  well  as  psychiatric  exclusions,  as  clarified  in  the  International  CFS  Working  Group  [60]. 

321 

322  Statistical  Analysis:  Brown’s  theoretical  approximation  [61]  of  Fisher's  statistics  was  used 

323  to  calculate  the  significance  of  alignment  between  experimental  data  and  a  given  model 

324  predicted  state  .  Fisher's  method,  a  meta-analysis  technique,  combines  probabilities  to 

325  obtain  the  overall  significance  of  a  set  of  P-values  obtained  from  independent  tests  of  the 

326  same  null  hypothesis.  The  combined  x2  statistic, 
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where  N  is  the  number  of  measureable  variables  and  p,  is  the  corresponding  P-values  under 
the  null  hypothesis,  has  a  x2  distribution  with  2N  degrees  of  freedom  assuming  that  the 
performed  tests  are  independent.  As  the  molecular  variables  of  the  endocrine  and  immune 
system  interact  with  one  another,  as  evidenced  by  the  above  connectivity  diagrams,  they  are 
not  independent.  As  a  result,  direct  application  of  this  test  statistic  is  invalid,  since  the 
assumption  of  independence  is  violated.  Brown  [61]  suggested  a  method  for  combining 
non-independent  tests.  If  the  tests  are  not  independent,  then  the  statistic  T0  has  mean  m  = 
2N  and  variance  (a2)  given  as, 

c>2=4N+2i=lN-l  j=i+lNcov(-2\npi,-2\npj)  (4) 

where  p,  and  py  are  the  P-values  for  each  test  and  the  covariance  (cov)  is  calculated  as, 

cov(-2\npi,-2\npj]=pij(3.25+0.75pij],  &  0<pij<lpij(3.27+0.71pij),  -0.5<pij<0 

(5) 

with  py  being  the  unadulterated  correlation  between  variable  i  and  variable  j.  Finally,  the 
overall  significance  P  of  a  set  of  non-independent  tests  is  calculated  using  the  statistic  T 
which  under  the  null  hypothesis  follows  the  central  x2  distribution,  where  T  =  T0/c  with  2N/c 
degrees  of  freedom  and  c  =  o2/4N. 

Here,  we  test  if  each  experimental  measure  aligns  with  a  given  model  predicted  state.  Our 
null  hypothesis  is  that  the  experimental  measures  do  not  align.  P-values  for  individual 
variables,  Pi,  are  calculated  using  two-sample  t-tests  between  ill  subjects  and  healthy 
controls.  Where  model  predictions  give  a  variable  as  high  (+1),  ‘right-handed’  one-tailed  test 
are  used,  whereas  a  ‘left-handed’  test  was  used  when  model  predictions  are  low  (-1),  to  give 
the  probability  of  obtaining  the  predicted  value  when  the  null  hypothesis  is  true.  For  the 
case  where  the  model  predicts  normal  behavior  for  a  variable  (0)  a  two-tailed  t-test  is  used. 
However,  the  p-value  from  the  two-tailed  test,  ptWo-taii,  gives  the  probability  that  there  is  an 
observable  difference  between  illness  and  control,  which  is  the  null  hypothesis.  To  rectify 

n 


354  this,  when  comparing  to  a  model  predicted  variable  of  0  we  take  the  P-value  to  be  Pi  =  1  - 

355  p two-taiu  giving  the  probability  of  obtaining  the  predicted  value  when  the  null  hypothesis  is  true. 

356 

357  All  cohort  data  was  normalized  using  a  Log2  transformation  before  T-tests  and  correlation 

358  calculations  were  performed.  The  unadulterated  correlation  values  py  between  two 

359  variables  i  and  j  were  calculated  in  healthy  subjects  as  the  pairwise  Pearson's  linear 

360  correlation  coefficient  between  variables  .  The  above-mentioned  experimental  data  was 

361  compared  against  model  predictions  based  on  the  five  measureable  variables,  namely 

362  TEST/EST,  CORT,  MR,  TICyt,  and  T2Cyt.  Where  model  variables  represent  an  aggregate 

363  set  of  markers  each  experimentally  measured  constituent  marker  was  compared  individually 

364  to  the  model  predicted  value.  For  example,  TICyt  is  composed  of  IL-2,  IFNy  and  TNF|3, 

365  therefore  3  individual  P-values  were  calculated  based  on  the  predicted  value  of  TICyt. 

366 

367  Results 

368  Stable  States  in  the  HPA  Models 

369  Application  of  the  discrete  state  representation  to  the  basic  stand-alone  HPA  model  (Figure 

370  1  A)  generated  27  system  states,  and  failed  to  produce  multiple  stable  states  (Figure  2). 

371  This  is  consistent  with  previous  ordinary  differential  equation  based  models  of  this  basic 

372  representation  of  the  HPA  axis  [21-26],  Discrete  state  representation  of  the  HPA-GR  model 

373  (Figure  1  B)  generated  243  system  states.  Of  these,  2  system  states  possessed  no 

374  outbound  edges  and  were  stable  attractor  steady  states  (Figure  2).  In  the  first  steady  state 

375  all  state  variables  assumed  nominal  values  whereas  the  second  steady  state  corresponded 

376  to  activation  of  state  variables  GRD  and  GR  with  suppression  of  ACTH  and  CORT.  Once 

377  again  this  solution  is  consistent  with  that  obtained  by  analysis  of  the  ordinary  differential 

378  equation  model  of  the  HPA-GR  system  proposed  by  Gupta  et  al.  [27]  and  Ben  Zvi  et  al.  [16], 

379 

380  Combining  the  HPA-GR  axis  with  the  HPG  axis  and  immune  system  (Figure  1  B-G) 

381  altogether  produced  4,782,969  system  states.  For  the  male  HPG  (model  a)  (Figure  1  C), 
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382  and  three  of  the  four  female  HPG  models  (models  b,  d  and  e)  (Figure  1  D,E,G)  five  steady 

383  states  were  identified  (Figure  2).  One  stable  state  is  characterized  by  nominal  values  for  all 

384  variables  (SSO),  which  corresponds  to  the  typically  normal  resting  state  of  the  system.  The 

385  first  alternate  state  (SSI )  displays  low  ACTH  with  high  GRD  and  GR,  while  the  second 

386  (SS2)  has  inhibited  innate  and  Thl  immune  responses  (low  ICell,  HR,  T 1  Cell,  and  TICyt), 

387  with  increased  Th2  activity  (high  T2Cell  and  T2Cyt).  The  third  stable  state  (SS3)  appears  to 

388  be  a  combination  of  SSI  and  SS2  with  low  ACTH,  ICell,  HR,  T 1  Cell  and  TICyt,  and  high 

389  GRD,  GR,  T2Cell  and  T2Cyt.  The  final  state  (SS4)  presents  with  hypercortisolism, 

390  suppressed  TEST  and  a  shift  towards  the  Thl  immune  reponse  (low  T2Cell,  T2Cyt,  GnRH, 

391  LH/FSH  and  TEST/EST,  and  high  CORT,  GRD,  GR,  TICyt  and  TICell).  The  persistently 

392  low  CORT  state  seen  in  the  previous  stand-alone  HPA  models  of  Gupta  et  al.  [16]  and  Ben 

393  Zvi  et  al.  [27],  was  not  recovered  here.  Instead,  CORT  was  expressed  at  a  nominal  or  high 

394  value  for  all  predicted  states.  SSI  most  closely  resembles  the  results  of  Gupta  et  al.  [27], 

395  and  Ben  Zvi  et  al.  [16],  however  these  previous  models  only  considered  a  single  regulator  of 

396  CORT,  namely  ACTH.  The  lack  of  a  predicted  hypocortisolic  state  in  SSI  here  can  be 

397  attributed  to  the  interplay  of  multiple  regulators  of  CORT  (ACTH,  HR,  TEST/EST,  and 

398  T 1  Cyt).  Inclusion  of  additional  regulators  is  not  expected  to  further  alter  this  state. 

399 

400  In  the  final  female  HPG  model  (model  c)  (Figure  1  F),  corresponding  to  the  ovulation  phase, 

401  these  same  five  states  were  recovered  along  with  six  new  additional  states  (Figure  2).  In  the 

402  first  three  additional  states  the  HPA  axis  and  innate  immune  response  are  suppressed  with 

403  low  CRH,  ACTH,  CORT,  ICell  and  HR,  while  the  HPG  and  anti-inflammatory  response  are 

404  raised  with  high  T2Cell,  T2Cyt,  GnRH,  LH/FSH  and  EST.  The  difference  between  the  three 

405  states  is  noted  in  the  level  of  glucocorticoid  receptor  response,  GR  and  GRD,  which  together 

406  take  values  of  low  (SS5),  nominal  (SS6)  and  high  (SS7).  The  remaining  three  additional 

407  states  all  give  suppressed  HPA  (CRH,  ACTH,  and  CORT)  and  lowered  TICell  activity,  with 

408  high  HPG  activity  (GnRH,  LH/FSH  and  EST),  and  are  again  differentiated  by  their 

409  glucocorticoid  receptor  levels  (GR,  GRD):  low  (SS8),  nominal  (SS9)  and  high  (SS10). 
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411  Overall,  inclusion  of  the  simplified  immune  system  and  the  HPG  works  to  regulate  CORT 

412  levels  in  the  HPA  axis.  The  male  HPG  (HPG  model  a),  and  the  majority  of  female  HPG 

413  configurations  (HPG  models  b,  d  and  e),  serve  to  produce  either  nominal  values  of  CORT, 

414  with  the  potential  of  a  shift  towards  Th2  activation  (SS2  and  SS3),  or  a  hypercortisolic  state 

415  with  low  TEST/EST  and  a  shift  towards  Thl  (SS4).  Only  connections  associated  with  the 

416  female  gender  (HPG  model  c)  were  responsible  for  the  emergence  of  a  natural 

417  hypocortisolic  state  (SS5  -  SS10).  This  hypocortisolic  state  comes  with  high  EST  and  may 

418  have  a  shift  towards  Th2  activation  in  the  immune  system. 

419 

420  Comparison  of  GWI  and  CFS  to  Predicted  States 

421  Application  of  Brown’s  meta-analysis  method  allowed  for  the  calculation  of  a  combined  P- 

422  value  comparing  the  experimental  data  with  the  predicted  stable  states,  allowing  for  the 

423  alignment  between  different  predicted  stable  states  to  be  ranked.  As  experimental 

424  measures  allowed  for  comparison  with  only  five  variables  (TEST/EST,  CORT,  HR,  TICyt, 

425  and  T2Cyt)  several  of  the  predicted  stable  states  resulted  in  the  same  experimental  profile 

426  and  resulting  combined  P-value  despite  being  distinct  states  (e.g.  SSO  and  SSI  both  show 

427  nominal  values  for  the  five  measureable  variables). 

428 

429  To  compare  to  our  model  the  difference  between  steroid  and  cytokine  levels  recorded  in 

430  male  Gulf  War  veterans  with  GWI  and  HCs  were  compared  to  the  steady  state  values 

431  predicted  by  the  male  variant  of  the  HPA-GR-lmmune-HPG  model  (model  a).  Comparison 

432  to  the  nominal  states  (SS0/SS1)  showed  poor  alignment,  Psso/ssi  =  0.82,  suggesting  that  the 

433  GWI  profile  cannot  be  considered  the  same  as  nominal  behavior.  Alignment  with  states 

434  presenting  a  shift  towards  Th2  immune  activation  (SS2/SS3)  showed  better  alignment, 

435  P ss2/ss3  =  0.38,  although  with  low  significance.  The  final  state,  displaying  hypercortisolism, 

436  low  TEST  and  a  shift  towards  Thl  immune  activation  (SS4),  yielded  the  best  alignment,  PSs4 

437  =  0.30,  again  however,  with  a  low  overall  significance. 
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439  The  difference  between  steroid  and  cytokine  levels  of  female  CFS  subjects  and  HCs  were 

440  compared  to  the  steady  state  values  predicted  by  the  female  variants  of  the  HPA-GR- 

441  Immune-HPG  models  (model  b-e).  Again,  alignment  with  states  presenting  nominal 

442  changes  in  measureable  variables  (SS0/SS1)  was  poor,  Psso/ssi  =  0.83,  supporting  that  CFS 

443  is  distinctly  different  from  normal  behavior.  The  Th2  shifted  immune  profile  states 

444  (SS2/SS3)  showed  a  significant  alignment,  Pss2/ss3  =  0.04,  suggesting  Th2  activation  in 

445  CFS.  This  is  further  supported  by  low  alignment  with  the  Thl  immune  activated  state,  with 

446  hypercortisolism,  and  low  EST  (SS4),  PSs4  =  0.28.  Improved  alignment  is  seen  in  states  with 

447  a  shift  towards  Th2,  coupled  with  hypocortisolism,  and  high  EST  (SS5/SS6/SS7),  Psss/ss6/ss7 

448  =  0.02,  suggesting  that  these  features  contribute  to  the  CFS  profile.  This  is  also  supported 

449  by  low  alignment  with  states  only  presenting  hypocortisolism  and  high  EST  with  no  immune 

450  activation  (SS8/SS9/SS10),  Pss8/ss9/ssio  =  0.60. 

451 

452  Discussion 

453  The  existence  of  multiple  stable  states  is  a  prime  characteristic  of  systems  incorporating 

454  feedforward  and  feedback  mechanisms,  and  plays  a  critical  part  in  guiding  the  complex 

455  dynamics  observed  in  biology.  These  alternate  stable  regulatory  regimes  occur  due  to  the 

456  feedforward  and  feedback  mechanisms  within  the  system  and  may  allow  escape  routes  for 

457  survival  of  an  insult  and  provide  support  in  the  medium  or  long-term  to  what  is  equivalent  to 

458  an  uneasy  cease-fire  or  adaptive  compromise.  An  example  of  such  compromises  in 

459  functional  status  in  exchange  for  survival  include  vasovagal  response  to  decreased  blood 

460  pressure  and  syncope  (“fainting”)  [62],  From  an  evolutionary  perspective  it  would  be 

461  advantageous  for  a  pathogen  to  establish  an  adaptive  relationship  with  the  host.  As  naturally 

462  occurring  alternate  states  of  homeostasis  are  inherently  stable  exploiting,  these  regimes 

463  could  be  an  advantageous  way  for  a  pathogen  to  establish  long-term  chronic  infection,  in 

464  essence  using  the  body’s  own  homeostatic  drive  to  maintain  the  status  quo.  To  explore  this 

465  hypothesis,  we  constructed  a  simple  but  integrated  model  incorporating  three  of  the  body’s 
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466  major  regulatory  axes:  the  HPA,  the  HPG  and  the  immune  system.  Modeling  the  dynamic 

467  properties  of  these  complex  systems  presents  a  significant  challenge,  as  much  of  the 

468  detailed  information  describing  in  vivo  kinetics  in  humans  is  unavailable.  However,  there  is  a 

469  very  significant  body  of  connectivity  data  describing  the  interactions  between  the  molecular 

470  and  cellular  elements  of  these  biological  systems.  To  make  use  of  this  wealth  of  information 

471  we  have  applied  a  discrete  state  representation  to  the  neuroendocrine  immune  system 

472  based  solely  on  the  biological  connectivity  found  in  the  literature  and  a  set  of  ternary  logical 

473  rules.  Using  a  discrete  logic  methodology  proposed  by  Thomas  [30],  we  demonstrated  that 

474  the  inclusion  of  feedforward/feedback  loops  leads  to  multiple  stable  states.  Indeed,  addition 

475  of  the  positive  feedback  loop  regulating  glucocorticoid  receptor  dimerization  (GR-GRD)  to  a 

476  basic  model  of  the  HPA  axis  generated  an  alternate  homeostatic  state  characterized  by  high 

477  receptor  expression  and  low  circulating  cortisol  levels,  a  result  found  previously  by  Gupta  et 

478  al.  [27]  and  Ben  Zvi  et  al.  [16]  using  differential  equation  based  models.  So  dependent  is  the 

479  natural  emergence  of  these  states  on  the  regulatory  wiring  that  inclusion  of  this  receptor 

480  dimerization  in  a  more  complex  HPA4mmune-HPG  models  resulted  in  the  disappearance  of 

481  this  alternate  hypocortisolic  state  through  compensatory  effects  of  these  axes.  Only  when  all 

482  three  interacting  axes  were  included  was  an  alternate  hypocortisolic  condition  recovered. 

483  Therefore  while  simple  models  require  the  inclusion  of  positive  receptor  feedback  dynamics 

484  to  produce  mutlistability,  these  effects  become  inherent  in  more  coarse,  but  comprehensive 

485  regulatory  circuits,  and  receptor-level  feedback  becomes  less  of  a  contributor  in  the  support 

486  of  multiple  attractor  states.  Coarse-grained  but  comprehensive  models  may  suffice 

487  therefore  in  capturing  physiologically  relevant  and  clinically  verifiable  response  dynamics. 

488 

489  Our  analysis  of  these  coarse  grained  models  spanning  across  multiple  regulatory  axes 

490  highlighted  the  important  role  of  gender  in  supporting  a  persistent  hypocortisolic  condition. 

491  Due  to  the  suppressive  actions  of  the  male  gonadal  system  in  regulating  itself  and  the  HPA 

492  axis,  a  low  cortisol  steady  state  is  never  available  to  the  male,  at  least  theoretically  at  this 

493  level  of  detail.  In  women  however,  the  combined  effect  of  EST  and  PROG  on  the  HPA  still 
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494  remains  somewhat  inconsistent  [34,63]  owing  to  the  varying  effects  of  these  hormones 

495  during  and  after  the  menstrual  cycle.  EST  is  generally  believed  to  stimulate  the  HPA  axis 

496  during  the  menstrual  cycle  [63-65],  however  evidence  indicates  that  in  perimenopausal, 

497  menopausal  or  ovariectomized  women  the  HPA  axis  response  is  inversely  correlated  with 

498  plasma  EST  levels  suggesting  an  inhibitory  effect  [65,66].  This  suggests  that  sex  hormone 

499  regulation  may  change  in  feedback  polarity  and  act  as  both  inhibitor  and  activator  of  the 

500  HPA  axis.  For  this  reason  HPA-HPG  interaction  in  women  will  in  theory  readily  support  the 

501  presence  of  a  stable  hypocortisolic  condition  when  HPG  axis  regulation  inhibits  the  HPA  axis 

502  while  stimulating  itself. 

503 

504  In  addition  to  sex  hormone  regulation,  interaction  with  the  immune  system  also  appears  to 

505  play  a  significant  role  in  determining  abnormal  cortisol  levels.  In  our  coarse-grain  models, 

506  cortisol  exerts  a  suppressive  action  on  the  innate  immune  system  and  the  Thl  adaptive 

507  immune  response.  Conversely,  positive  feedback  by  certain  components  of  the  immune 

508  system  promotes  increases  in  cortisol  levels,  which  support  a  hypercortisolic  steady  state. 

509  While,  inclusion  of  the  glucocorticoid  receptor  dimerization  (GR-GRD)  in  these  models 

510  yielded  additional  steady  states,  it  did  not  result  in  any  significant  changes  to  the  profile  in 

511  regards  to  cortisol  levels.  Combining  the  actions  of  HPA,  HPG  and  immune  regulation 

512  supported  the  existence  of  a  stable  hypercortisolic  state  in  all  models  of  men  and  women 

513  while  a  persistent  hypocortisolic  state  was  available  only  in  women  and  only  under  certain 

5 14  modes  of  HPG  regulation.  Once  again,  while  the  inclusion  of  the  GR-GRD  receptor 

515  dimerization  in  this  overarching  model  yielded  additional  steady  states,  it  did  not  result  in  any 

516  significant  changes  to  the  homeostatic  profiles. 

517 

518  These  findings  suggest  that  abnormally  high  levels  of  cortisol  and  adaptive  immune 

519  activation,  in  this  case  Thl ,  may  be  perpetuated  under  certain  conditions  by  the  system’s 

520  own  homeostatic  drive.  This  prediction  of  persistent  and  stable  Thl  activation  is  consistent 

521  with  evidence  of  anomalies  in  immune  signaling  in  GWI  [55,67,68],  Skowera  et  al.  measured 
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522  intracellular  production  of  cytokines  in  peripheral  blood  and  found  ongoing  Thl-type  immune 

523  activation  in  symptomatic  Gulf  War  Veterans  compared  to  healthy  counterparts  [67],  More 

524  recent  work  confirmed  this  finding  while  also  suggesting  that  this  may  occur  in  the  more 

525  complex  context  of  a  mixed  Thl  :Th2  response  [55],  something  not  captured  by  the  simple 

526  immune  model  used  here.  Though  we  were  unable  to  find  documented  reports  of  lower 

527  testosterone  levels  in  GWI  beyond  the  experimental  data  presented  here,  a  large  study  of 

528  gulf  war  veterans  in  the  UK  found  increased  risk  of  fertility  problems  in  this  population  [69], 

529  suggesting  a  possible  relation. 

530 

531  In  much  the  same  way,  conditions  involving  hypocortisolism  and  a  Th2  shift  may  also  be 

532  perpetuated  at  least  in  part  by  the  natural  homeostatic  regulatory  programming.  In  this  case 

533  the  homeostatic  program  may  be  driven  by  sex  steroid  suppression  of  the  HPA  axis  and 

534  promotion  of  HPG  function  coupled  with  the  mutual  inhibition  between  the  Thl  response  and 

535  function  of  the  gonadal  axis,  a  configuration  seemingly  available  only  to  female  subjects  in 

536  our  models.  This  would  suggest  that  the  hypocortisolism  seen  in  diseases,  such  as  CFS  [70- 

537  72],  could  be  a  result  of  the  complexity  afforded  by  the  interaction  between  the  HPA, 

538  immune  and  HPG  axes  in  female  subjects.  Indeed  model  predictions  describing  such  an 

539  alternate  homeostatic  state  in  women  aligned  with  our  experimental  results  from  CFS 

540  subjects,  and  is  consistent  with  previous  findings  of  Th2  activation  in  CFS  (Brenu  et  al. , 

541  201 1 ,  Nakamura  et  al.,  2010  and  Natelson  et  al.,  2005,  Broderick  et  al.,  2010).  This 

542  alignment  with  a  naturally  occurring  homeostatic  conditions  may  explain,  at  least  in  part,  the 

543  biased  prevalence  of  such  persistent  diseases  in  women  [73-78],  Indeed,  these  authors 

544  report  that  approximately  70%  of  observed  CFS  patients  are  women.  Additionally,  the 

545  prevalence  of  CFS  in  the  40-49-year-old  age  range  [78],  and  the  higher  prevalence  of 

546  gynecological  conditions  and  gynecological  surgeries  in  women  with  CFS  [79]  supports  the 

547  evidence  that  HPA  suppression  by  estradiol  appears  more  likely  in  perimenopausal, 

548  menopausal  or  ovariectomized  women  [65,66],  Interestingly,  as  many  as  1  in  3  CFS 

549  subjects  have  reported  symptom  relief  during  pregnancy  [80],  The  normal  trend  in 
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pregnancy  towards  increased  cortisol  levels,  especially  in  the  third  trimester,  might  be  a 
contributing  factor  that  would  support  the  key  involvement  of  sex  hormone  regulation 
proposed  by  our  analysis  [81].  While,  in  normal  pregnancy  this  increase  in  cortisol  typically 
coincides  with  an  increase  in  cortisol-binding  globulin  (CBG)  maintaining  the  level  of  free 
cortisol,  CBG  genetic  variants  in  CFS  have  the  potential  to  alter  normal  CBG  function 
[82,83], 

While  certainly  more  comprehensive  than  their  predecessors,  these  models  remain  relatively 
coarse  representations  of  the  interplay  between  the  endocrine  and  immune  systems. 

This  is  particularly  true  of  immune  model  granularity,  especially  when  one  considers  the 
complex  signaling  network  supported  by  immune  cells  as  well  as  other  immune-sensitive 
cells  [84],  The  important  role  of  key  neurotransmitters  linking  the  central  nervous  system 
with  the  HPA  axis  and  the  immune  system  was  also  under-represented  in  this  first 
generation  of  models.  For  example,  norepinephrine  and  epinephrine  stimulate  the  (32- 
adrenoreceptor-cAMP-protein  kinase  A  pathway  inhibiting  the  production  of 
Thl/proinflammatory  cytokines  and  stimulating  the  production  of  Th2/anti-inflammatory 
cytokines  causing  a  selective  shift  from  cellular  to  humoral  immunity  [85,86],  Additionally, 
lymphocytes  express  most  of  the  cholinergic  components  found  in  the  nervous  system. 
Lymphocytes  may  be  stimulated  by,  or  release,  acetylcholine  thus  constituting  an  immune 
regulating  cholinergic  system  secondary  to  the  nervous  system  [87],  Another 
neurotransmitter,  neuropeptide  Y  (NPY),  also  serves  as  a  powerful  immune  modulator  [88] 
and  has  recently  been  shown  to  play  a  role  in  CFS  [89].  These  components  are  without 
question  important,  however  based  on  our  initial  observations  from  this  piecewise  analysis 
we  expect  that  increased  detail  will  lead  to  the  emergence  of  additional  response  programs 
rather  than  the  elimination  of  attractors  found  here. 

As  these  models  are  based  on  currently  documented  knowledge  of  human  physiology  and 
regulatory  biochemistry  they  are  necessarily  incomplete.  Nonetheless  the  simple  models 
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578  presented  here  illustrate  the  importance  of  an  integrative  approach  to  understanding 

579  complex  illnesses.  Further  refinement  of  the  model  to  include  more  detailed  description  of 

580  interactions  within  and  between  the  HPA,  HPG  and  immune  systems  could  extend  its 

581  applicability  to  other  illnesses  as  would  the  incorporation  of  other  key  systems  such  as  the 

582  brain  and  central  nervous  systems.  Yet,  even  with  the  coarse-grained  co-regulation 

583  networks  investigated  we  found  numerous  stable  resting  states  that  differ  significantly  from 

584  normal  and  were  indicative  of  complex  and  persistent  regulatory  imbalances.  Findings  such 

585  as  this  support  the  use  of  an  alternate  model  for  disease,  one  which  is  not  necessarily 

586  associated  with  failure  of  individual  components,  but  rather  with  a  shift  in  their  coordinated 

587  actions  away  from  normal  regulatory  behavior.  Response  to  exercise  and  other  stressors 

588  has  the  potential  to  be  very  different  in  these  new  regulatory  regimes.  This  is  something  that 

589  we  have  observed  firsthand  in  our  work  with  human  GWI  and  CFS  subjects  [90], 

590 

591  Finally,  when  considering  alignment  with  the  experimental  data  presented  here  for  CFS  and 

592  GWI,  it  is  important  to  remember  that  it  was  never  our  hypothesis  that  these  illnesses 

593  resulted  solely  from  the  actions  of  homeostatic  drive.  Instead  we  proposed  that  homeostatic 

594  drive  might  be  a  significant  contributor  to  the  persistence  of  illness  mechanisms.  Because 

595  these  naturally  occurring  regimes,  once  instantiated,  provide  an  alternate  stable 

596  homeostasis  resistant  to  change,  it  may  offer  fertile  ground  in  support  of  many  chronic 

597  pathological  processes.  The  alignment  of  several  immune  and  endocrine  markers  modeled 

598  here  with  experimental  data  from  CFS  and  GWI,  two  chronic  conditions,  would  support  at 

599  least  partial  involvement  of  the  body’s  own  homeostatic  drive  in  facilitating  the  perpetuation 

600  of  these  conditions.  This  may  promote  resistance  to  therapy  and  the  natural  regulatory 

601  barrier  to  change,  even  positive  change,  should  at  least  be  considered  in  the  design  of 

602  robust  treatment  avenues. 
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Tables 


Table  1:  Ternary  HIGH/LOW  PASS  operator 


A  VB 

B  =  -1 

B  =  0 

B  =  1 

A  =  -1 

0 

0 

-1 

A  =  0 

0 

0 

-1 

A  =  1 

1 

1 
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Table  2:  Ternary  OR  opera 

tor 

AvB 

B  =  -1 

B  =  0 

B  =  1 

A  =  -1 

-1 

0 

1 

A  =  0 

0 

0 

1 

A  =  1 

1 

1 
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Table  3:  Ternary  NOT  operator 


A 

-A 

-1 

1 

0 
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1 

-1 

Figure  Legends 

Figure  1:  Standard  and  extended  HPA  models.  (A)  Standard  HPA  model.  (B)  HPA-GR 
model  of  Gupta  et  al.  [27],  Integrated  models  (C)  HPA-GR4mmune-HPGa  for  males,  and  (D) 
HPA-GR4mmune-HPGb,  (E)  HPA-GR4mmune-HPGc,  (F)  HPA-GR4mmune-HPGd,  and  (G) 
HPA-GR4mmune-HPGe  for  females. 


Figure  2:  Steady  states  of  standard  and  extended  HPA  models.  White  -  nominal  state  (0); 
Green  -  high  state  (1 );  Red  -  low  state  (-1 );  Grey  -  N/A  to  the  model. 
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Figure  2 
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Appendix  B: 

Updated  scope  of  work  (SOW)  submitted  as  part  of  request  for  transition  of  award  to  Nova  Southeastern 
University. 
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Revised  Statement  of  Work  (SOW)  -  Relocation  to  Nova  Southeastern  University  (NSU),  FL 


The  following  SOW  has  been  updated  to  reflect  current  status  of  project  consistent  with  the  annual  progress  report  submitted  in 
September.  2012.  Described  are  the  remaining  activities  required  for  completion.  These  will  now  be  conducted  by  the  Broderick 
group  from  its  new  home  institution:  Nova  Southeastern  University,  Fort  Lauderdale,  Florida. 

In  support  of  Dr.  Broderick’s  transition  Nova  Southeastern  University  is  entering  into  a  service  agreement  with  the  Center  for 
Computational  Sciences  (CCS)  at  the  University  of  Miami,  which  will  serve  as  the  principal  high-performance  computing  resource 
for  the  Broderick  group  from  this  point  forward.  The  latter  will  continue  to  use  the  University  of  Alberta’s  WestGrid  high- 
performance  computing  platform  during  the  transition  period  in  order  to  ensure  continuity  of  the  work. 

Taskl.  Evaluate  and  select  agent-based  simulation  environment.  Completed. 

Task  2.  Define  and  encode  immune  cell  populations  and  interaction  rules.  Completed. 

Task  3.  Refine  HPA  axis  model  and  integrate  with  immune  model.  Completed. 

MILESTONE  I:  Completion  and  release  of  validated  model  combining  ODE  representation  of  the  HPA  axis  and  a  discrete 
population-based  model  of  the  immune  system.  Target  date:  Completed. 

Extensions  to  original  Task  2,  3. 

3. a.  Extension  of  circuit  model  of  neuro-inflammatory  cascades.  In  an  extension  of  the  original  mandate  for  Task  2,  the 

circuit  logic  approach  is  also  being  applied  to  model  inflammatory  processes  occurring  in  the  brain  and  involving  the  cell 
types  and  immune  signaling  specific  to  this  physiological  compartment. 

3.b.  Extension  to  sex  hormone  and  thyroid  axes.  Task  3  has  also  been  extended  beyond  the  original  mandated  scope  to  now 
include  the  endocrine  axis  regulating  sex  hormones.  We  expect  to  integrate  thyroid  function  in  this  regulatory  circuitry  as 
well. 

Timeline  extended  mandate:  Months  1-4,  Year  3  (now  December,  2013) 

Site(s):  Nova  Southeastern  University 

Task  4.  Design  and  conduct  formal  sensitivity  and  multi-stability  analyses.  Completed 
Task  5.  Network  analysis  of  alternate  homeostatic  states.  Completed 

Extensions  to  original  Task  4,  5.  These  steps  will  be  repeated  in  the  analysis  of  the  extensions  to  the  model  proposed  in  3(a)  and 
3(b). 

Timeline  extended  mandate:  Months  1-4,  Year  3  (now  December,  2013) 

Site(s):  Nova  Southeastern  University 

MILESTONE  II:  Verification  of  hypothesis  that  GWI  symptoms  persist  because  the  endocrine-immune  system  now  occupies  and 
alternate  homeostatic  stable  point  and  engages  a  new  sub-optimal  stress  response  control  program. 

Target  date:  Month  4,  Year  3  (now  December,  2013);  Currently  70%  complete. 

Task  6.  Identify  and  deploy  large-scale  optimization.  This  involves  the  selection  of  the  best  algorithm  for  exhaustive  search  of 
intervention  possibilities.  We  expect  the  combined  endocrine-immune  system  to  present  multiple  stable  points  and  the 
landscape  describing  its  dynamic  response  to  be  complex.  As  a  result  standard  techniques  for  optimization  of  treatment  time 
course  would  terminate  their  search  in  the  first  region  where  treatment  performance  ceases  to  improve.  Overall  such  a 
treatment  may  be  quite  remote  from  that  available  in  the  neighboring  response  “valley”. 

Timeline:  Months  2-4,  Year  3  (now  October  -  December  2013) 

6. a.  Review  global  search  algorithms.  Review  latest  developments  in  evolutionary  programming  techniques  as  well  as 
hybrid  techniques  to  determine  the  most  suitable  search  algorithm.  Acquire  or  develop  code  and  deploy  on  CCS 
platform  and  test  on  logic  model  developed  in  4b. 

Timeline:  Month  1-2,  Year  3  (now  September  -  October  2013) 

Site(s):  Nova  Southeastern  University 

6.b.  Configure  simulation-based  optimization  scheme.  Configure  an  interface  that  evaluates  the  fitness  of  candidate 
interventions  by  repeatedly  launching  short  logic  model  simulations  as  it  conducts  its  search  for  the  most  robust 
treatment  course.  Test  and  deploy. 

Timeline:  Months  2-6,  Year  3  (now  October  -  February,  2014) 
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Site(s):  Nova  Southeastern  University 


Task  7.  Identify  candidate  treatment  courses  for  GWI.  Using  the  optimization  scheme  developed  and  deployed  in  Task  6,  launch 

optimization  runs  from  multiple  initial  conditions  of  endocrine-immune  status.  Assess  these  options  and  report. 

Timeline:  Months  5-12  Year  3 

7. a.  Define  and  encode  solution  fitness  criteria.  Identify  and  describe  mathematically  the  immune  and  endocrine  descriptors 
that  can  be  safely  changed  and  over  what  range  they  may  be  changed.  Incorporate  these  constraints  with  treatment  goals 
and  define  optimization  problem  formally. 

Timeline:  Month  3-5,  Year  3  (now  November,  2013  -  January,  2014) 

Site(s):  Nova  Southeastern  University 

7.b.  Search  for  broadly  applicable  candidate  treatment  courses.  Identify  a  set  of  initial  conditions  of  cytokine,  hormone  and 
immune  cell  abundance  and  launch  repeated  searches  for  optimal  treatments  from  these  points. 

Timeline:  Months  6-9,  Year  3  (now  February  -  May,  2014) 

Site(s):  Nova  Southeastern  University 

7.c.  Critically  assess  candidate  treatments.  Review  candidate  treatment  courses  and  assess  these  critically  based  on  efficacy 
and  minimal  invasiveness.  Propose  design  of  pilot  clinical  trials  for  evaluation  of  the  best  candidates. 

Timeline:  Months  10-11,  Year  3  (now  June  -  July,  2014) 

Site(s):  Nova  Southeastern  University 

7.d.  End  of  project  review  and  report.  Timeline:  Month  12,  year  3  (now  August,  2014). 
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